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SUMMARY 


The recent interest in radiofrequency (rf) induction heating as a 
means of producing a plasma has prompted work on theoretical techniques 
that can adequately model the physical processes occurring in the plasma. 
Efforts to solve the rf induction plasma problem have been hampered by a 
lack of complete input data reported by the experimenters, most notably 
the work coil geometry and the magnetic field intensity at the plasma 
radius . 

The numerical techniques and approximate analytical closed form 
solutions of the energy balance equation and Maxwell's equations used to 
describe the plasma require a knowledge of the electric and magnetic field 
values at the plasma radius and the coil geometry. The only comparisons 
of theoretical solutions with experimental results to date have been in- 
direct comparisons based on reported power requirements. 

Measurements of the temperature profiles in an rf argon plasma 
were made over magnetic field intensities ranging from 20 amp turns/cm to 
80 amp turns/cm. The results were compared with a one-dimensional numeri- 
cal treatment of the governing equations and with an approximate closed 
form analytical solution that neglected radiation losses. 

The average measured temperatures in the plasma compared well with 
the numerical treatment, though the experimental profile showed less of an 
off center temperature peak than predicted by theory. This may be a re- 
sult of the complex turbulent flow pattern present in the experimental 
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torch and not modeled in the numerical treatment. The radiation term 
cannot be neglected for argon at the power levels investigated. The 
closed form analytical approximation that neglected radiation led to tem- 
perature predictions on the order of 1000°K to 2000°K higher than measured 
or predicted by the numerical treatment which considered radiation losses . 

The governing equations are highly coupled to the electrical con- 
ductivity of the working fluid. An aerosol of submicron tungsten particles 
was added to the argon to effect a significant increase in the electron 
number density and temperature measurements were taken and compared with 
the theory. The temperatures were of the order of 1000 °K lower than for 

the cases of pure argon, and compared well with the theory at tungsten 

-3 -3 

partial pressures of 4.3 x 10 atm and 6,1 x 10 atm. A tungsten par- 

-2 

tial pressure of 3.6 x 10 atm led to a divergence between the theory 
and the experiment as the tungsten atom number density reached the point 
that electron-tungsten collisions could no longer be neglected in the 
theoretical treatment of the electrical conductivity. 

A parametric study was conducted to determine the effect of varia- 
tions in the inputs to the analysis which can be varied experimentally, 
and to determine the effect that the uncertainty in the transport proper- 
ties of argon has on the predicted temperature profiles. 
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CHAPTER I 
INTRODUCTION 

The recent interest in radiofrequency (rf) induction heating as a 
means of producing high density plasmas has prompted work on theoretical 
techniques that can adequately model the physical processes occurring in 
the plasma. The theoretical treatment must include the heat transfer 
problem and the determination of the electromagnetic field intensities 
inside the plasma. Previous theoretical work on modeling the d.c. arc 
plasma forms a basis for treating the energy balance problem, but these 
equations must now be treated simultaneously with the electromagnetic 
field equations as they apply to the induction arc. 

Efforts to solve the rf induction plasma problem have been hampered 
by a lack of complete input data reported by the experimenters. The work 
coil geometry and the magnetic field intensity at the plasma radius are 
needed as starting points for quantitative comparisons of spectroscopi- 
cally measured temperature profiles. Experimentally measuring the mag- 
netic field intensity at the plasma radius is difficult since the coil 
currents are high (of the order of 50-150 amperes, rms) and the frequen- 
cies are high (typically between 2 and 30 MHz). 

The development of the induction arc stems from work by Hittorf 1 

2 

on the low pressure electrodeless ring discharge. Babat is credited 
with inventing the induction arc in 1942 when he used a powerful vacuum 
tube oscillator to raise the pressure of a ring discharge up to atmospheric 
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level. T. B. Reed first demonstrated in 1961 that an "induction torch" 
at one atmosphere could be produced. He used a 26 mm outside diameter 
open ended quartz tube surrounded by a 3/16 inch copper work coil to pro- 
duce an argon plasma. Spectroscopic temperature profiles of argon plas- 

4 5 

mas were determined by Goldfarb and Dresvin, Schultz and Anderson, and 
Rovinskii, Gruzdev, et al #J ^ and the results of their work are presented 
in Figure 1, taken from Eckert 1 survey of the state-of-the-art of in- 

g 

duction plasmas. Stokes reported temperature profiles of an argon plasma 

at one atmosphere along with the power dissipated per unit length of work 
9 10 

coil which Eckert ’ used as a basis for comparison with the numerical 

solutions of Pridmore-Brown^ and his own closed- form analytical predic- 
9 12 

tions. ’ Stokes used a water-cooled quartz tube with a 2.3 cm inner 
radius and excited the gas at 3.8 MHz. The power dissipation was deter- 
mined from the measured heat flux into the cooling water and from the 
total radiation emitted. The applied magnetic flux was not reported. A 
wide range of magnetic flux intensities was used in the calculations 
until the predicted radiation and conduction values matched those re- 
ported in the experiment. The predicted temperature profiles were in 
good agreement with spectroscopically measured values though the experi- 
mental profiles were flatter, as seen in Figure 2. The experimental tem- 
perature profiles were determined from continuum intensity measurements 
at 4288 A. 

This dissertation describes work aimed at providing spectroscopi- 
cally determined radial temperature profiles in an rf produced argon 
plasma over a wide range of experimentally measured applied magnetic flux 



Temperat 



Radius, R, (cm) 


Figure 1. Measured Radial Temperature Distributions in Several 
Induction Plasmas of Different Radii 
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Figure 2 


Comparison Between Spectroscopically Determined Temperatures 
Numerical Calculations, and Closed Form Analytical Solutions 



5 


intensities and to compare the results with a representative numerical 
analysis based on the one-dimensional Elenbaas-Heller equation and with 
an approximate closed form analytical solution. A better judgment can 
then be made on the accuracy of these theoretical treatments. Additionally, 
the argon plasma is seeded with submicron tungsten particles to effect a 
significant change in the electrical conductivity of the gas and a com- 
parison is made between the experimentally observed changes in the tem- 
perature profiles with those predicted by existing theory. 
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CHAPTER II 
THEORETICAL STUDIES 
Introduction 

The heat balance at each position in an induction plasma arc is 

given by the Elenbaas-Heller equation which balances the heat generated 

through current dissipation with conduction, convection, and radiation 
13 14 

losses. * The heat balance equation is coupled to Maxwell's equations 

governing the electric and magnetic fields in the induction plasma and 

both sets of equations must be satisfied. 

Analytical and numerical attempts to solve the problem are hindered 

by a lack of experimental data that include all input conditions required 

for quantitative comparisons, most notably the magnetic field strength 

inside the work coil. Numerical methods are cumbersome and require a 

digital computer. Analytical studies have failed to produce quantitative 

predictions universally applicable to the variety of experimental condi- 

11 12 

tions and gases used in rf plasma work. Eckert and Brown * compared 

Q 

their numerical method to the experimental work of Stokes, but even in 
this case Stokes did not measure the ampere turns of the induction coil. 
This measurement is needed to determine the magnetic field strength inside 
the work coil. 

The transport properties of argon are examined in this chapter 
because the quantitative results of the solution to the plasma problem 
are dependent on the values chosen for the transport properties and the 



7 


theory serves as a basis for predicting the properties of a tungsten 
seeded argon plasma. 


Governing Equations 

Consider the plasma to be an infinite axially symmetric cylinder 
with an applied alternating axial magnetic field 


H z = H o exp(i2m>t) 


( 1 ) 


where v is the frequency of the alternating current I in the work coil as 
shown in Figure 3. 

Maxwell's Equations 

The relevant Maxwell equations'^ are 


V 


X E = 


6B 6H 

& t at 


( 2 ) 


V X H = J + 




bt 


6t 


(3) 


J = aE 


( 4 ) 


where E and H are the induced electric and magnetic field intensities, B 
is the magnetic induction, D is the electric displacement current, J is 
the electric current density, and a is the electrical conductivity of the 
plasma. Equation (4) assumes a purely resistive plasma, i.e., there is 
no imaginary component of the conductivity. In the absence of polariza- 
tion or magnetization. 
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B = (jiH 

(5) 

D = sE 

(6) 


where |j, is the permeability and e the permittivity of the gas. The 
Poynting vector S, defined by S = E X H, represents the flux of electro- 
magnetic energy and points in the direction of propagation of the wave. 

1.6 

Following the treatment of Brown, et al. as modified by Freeman 
and Chase, ^ equation (2) expressed in cylindrical coordinates for an 
electric field with component E only and magnetic field with axial com- 

V 

ponent H only becomes: 

Z 


dE e 

E e + r 3F 



(7) 


or, from equation (1) 


r to rE 9 = • i2Tr ^ K z (8) 

Similarly, equations (3) and (4), in the absence of displacement 
currents, can be expressed in cylindrical coordinates for the applicable 
electromagnetic field components 


dH 

dr” J e = ' aE e 


(9) 


The assumption of a good conductor (i.e,, no displacement currents) 


is equivalent to requiring 
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a » e2nv (10) 

This is a valid assumption for argon at 4 MHz over the temperature 
range investigated as can be seen from the electrical conductivity values 
for argon presented later. 

Equations (8) and (9) may be combined to eliminate either H or E.. 

Z U 

Differentiating equation (8) and substituting into equation (9), the equa- 
tion for the induced electric field is obtained 

d2E e i dE / n 

— + 7 d ” ‘ i i2TTV ^ c + ~2J E © = 0 (u) 

dr r 

If E a is divided into its real and imaginary parts, equation (11) results 
in two equations 


d 2 E dE E 

71 “ + 7 d7 2 + E I " 0 

dr r 


d 2 E dE E 

~T + 7 iT - -i - e r * 0 

dr r 


(12) 


(13) 


The rate at which electromagnetic energy in a unit volume is dis- 
sipated can be arrived at from equations (2) and (3). Scalarly multiply- 
ing equation (2) by H and equation (3) by -E and adding the two equations 
produces 

H*(V x E) - E-(V X H) = - ^ (txH 2 + eE 2 ) - J-E 


(14) 



11 


Now using the vector identity 


7.(E X H) = H- (V x E) - E-(V x H) 


(15) 


to get 

V.(E X H) * - j^r CiH 2 + eE 2 ) - J.E 


(16) 


and integrating over a unit volume dr, and applying Green's theorem to 
the left hand side of equation (16), one obtains 


f (E X H)-ds = - -^ [ (e£ 2 + p,H 2 )dr - J J.EdT (17) 

s 0t *V 

The first term on the right hand side represents the time rate of 
decrease of stored electromagnetic energy in the volume and the second 
term on the right hand side represents losses by Joule heating or gains 
by an applied electric field. The left hand side then must represent the 
net outward flow of energy across surface s and leads to the definition 
of the Poynting vector S 


S = E X H (18) 

which represents the flux of electromagnetic energy. For steady state, 
Poynting’ s theorem is then 

V.S = - J.E 


(19) 
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or from equation (4) and for our geometry 

V.S = - o|e| 2 (20) 

If E is divided into its real and imaginary parts, E^ and E^, the 
time averaged root mean square value of the electric field in the 0 direc- 
tion is then 


B e . 4F CEl > + ^ 


( 21 ) 


and the power input per unit volume for use in the steady state energy 

2 

balance equation is then aE Q or 


aE 0 2 = i a (E x 2 + E^) 


( 22 ) 


Energy Balance Equation 

The energy balance equation equates the time rate of change of 
energy in a unit volume with the divergence of the heat flux, the power 
radiated, and the electromagnetic power dissipated by the applied electric 
field 

p C p dI = V ’ XVT + aE e 2 - * (23) 

DT 

where ^ stands for the total time derivative of the temperature, T, p is 
the plasma density, C p the specific heat at constant pressure, \ is the 
thermal conductivity, and § is the radiation source strength. The second 
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term on the right side of the equation represents the dissipative compo- 
nent of the electromagnetic energy. 

Expanding equation (23) in cylindrical coordinates and assuming a 
constant pressure axisymmetric plasma, one obtains 


+ £[(© l+ (©W-*-o 


(24) 


where and represent the radial and axial gas velocities. 

Assuming steady state conditions and rearranging the terms, one 
obtains the two-dimensional energy balance equation 


1 (SlI . 1 bT , . dX 17 bT\ 2 /6T\ 2 “I 

V 2 r br ^ " Lfey + fey J 


(25) 


- pC (u + u + crE. 

K p \ r 5r z bz/ 6 


- § = 0 


The energy balance equation may be reduced to one-dimensional form 
if it is assumed that the temperature gradients in the axial direction are 
very small compared to the temperature gradients in the radial direction 



I ±_ 

r dr 


r 



H T 4 

pC u “ + cE ' 
K p r dr 6 


$ = 0 


(26) 
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The radial flow velocity, u^, is assumed zero in the one-dimensional 

9 18 19 20 

analytical treatments of Eckert * and Rovenskii, et al., 9 and in the 

11 21 

one-dimensional numerical treatments by Brown, and Hiller and Ayen. 

22 23 

The work of Chase, ’ however, shows that the rf field imposes a magnetic 
pressure gradient that causes the charged particles to preferentially mi- 
grate toward the center of the plasma, and through collisions give a net 

inward radial velocity that forms a double vortex flow. Experimental work 
22 

by Chase on an rf plasma with a configuration similar to the one used in 
this research shows the radial velocity to be of the order of magnitude of 
10% of the net axial flow velocity. 

Equation (26) may be linearized by introducing the heat conduction 
potential defined by 

T 

s « J X dT (27) 

o 


With the aid of equation (27) , equation (26) becomes 


d 2 S 1 dS 

. 2 + r dr 
dr 


° C p U r dS 2 

X dr + E e 


$ = 0 


(28) 


If E is divided into its real and imaginary parts as shown in 
0 

equation (21) , equation (28) becomes 


d 2 S 1 dS 

, 2 + r dr 
dr 


pC u 
M P r 


< v > - 4 - 0 


( 29 ) 


and equations (12,13, and 29) along with the continuity equation in 
cylindrical coordinates 
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— 4 — (rpu ) = 0 

r dr H r' 


( 30 ) 


and boundary conditions: 


r = 0 


wall 



s = 

J * dT 

dr 

dE R 

O 

Er = 0 

r = R ^ 
dr 

o 

II 

+ 

W 

M 

II 

O 

dE T 

E t 


I 

dr 

+ -7 = - 2ttmvH(R) 


(31) 


form a one-dimensional system of equations that describe the plasma. 

Solution Techniques 

Eckert ^ provides an up to date review of the reported techniques 
used for numerical solutions to the equations governing the induction arc 
as well as the assumptions necessary for approximate closed form solutions. 
He pointed out that lack of complete input data from reports of experimen- 
tal work makes quantitative comparisons difficult and that experimenters 
did not attempt to minimize two-dimensional or convection effects in their 
work. 

24 

Keefer at the University of Florida, developed a one-dimensional 

25 

numerical technique programmed under the direction of Sprouse at the 
Arnold Engineering Development Center. This program (called SUPERSIX) was 
obtained and converted to run on the Univac 1108 at Georgia Tech to permit 
a comparison of the experimental work reported in this dissertation with 
a numerical analysis of the governing equations. 

The analysis reduces the sixth order system of equations given by 
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equations (12, 13, 29, 30) to six first order differential equations as 
follows* Defining 



o 


(32) 


where is the input wall temperature and 


y 2 



(33) 


y, = 


y 4 



y s = E t 


y 6 



where x is the non-dimensional radial position x = r/R, one obtains from 
equations (12, 13, 29, 30) 

dy i 

3T - y 2 < 34 > 

dy 2 /*P u a C p< y l> 1\ r2 ' ,( 5'i ) , 2 2, R2{ <yi> 

dT “ l Tx ■ x) y 2 - “3— < y 3 + y 3 > + — 

w w 


(continued) 
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^3 = 1 

dx x y 4 

dy 4 1 2 

i = - y - 2R xqj.2TTvy, 
dx x 3 5 



X y 6 


dy 6 1 2 
— - - i y + 2R xou,2TTvy- 
dx x J 5 -i 


where R is the plasma radius with boundary conditions 

x = 1 : y x = 1; y 5 + y 6 = - 2 ttj,vRH(R) ; y 3 + - 0 (35) 

x = o : y 2 = y 3 = y 5 ~ ° 

A Runge-Kutta integration scheme is used to solve the system of equations 
(34,35) with the required inputs being frequency, wall temperature, plasma 
radius, and magnetic field strength at the outer surface of the plasma. 

The radial flow velocity at the wall is input through the term pu r R. Pro- 
visions are made in the program for parametric analysis using the last 
plasma profile as the initial guess for a new case where only one input 
parameter is slightly varied. The analysis does not take axial convection 
into account, as can be seen from the governing equations. Typical run 
times on the Univac 1108 for convergence to the first case are 20-30 
minutes • 

19 20 

The analytical work of Gruzdev, Rovenskii, and Sobolev ’ provides 
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a simplified closed form solution to the one-dimensional problem that is 
useful for comparison with the numerical method and the results of this 
research. By considering E and H as root mean square values such that the 
value of the Poynting vector becomes E-H and neglecting radial convection 
and radiation, the one-dimensional energy balance equation for steady 
state becomes 


V-XVT + 7 . EH =* 0 


or in cylindrical coordinates 


il( 

r dr \ 


rEH - r X ~) = 0 


Since both fluxes vanish at r = 0 


dT 

dr 


= EH 


(36) 


(37) 


(38) 


From equation (9) 


dT _ H dH 

dr a dr 


(39) 


Integrating equation (39) from the outer boundary r = R where the temper- 
ature equals the wall temperature, T^, and the magnetic field has a value 
at the boundary of H„ 
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If the magnetic field as a function of r is known, the temperature at 

any point may be calculated* If H(0) « which is a good assumption 

26 

according to the experimental work of Eckert and the theoretical predic- 
tions of the numerical analysis, then the maximum temperature (which occurs 
at the center of the plasma in the case of no radiation) may be approxi- 
mated from the magnetic field at the plasma radius from 


C p T(0) -A 

h r = L 2 J a(T) X(T) dT J (41) 

T w 

The choice of T does not affect the final T(0) for a given H for 
w k 

argon as long as T w < 4000 °K because of the very low values of a and \ 
below about 5000°K. 

Transport Properties 

The values chosen for the transport properties of argon, namely 
the electrical and thermal conductivities and radiation source strength, 
bear directly on the quantitative results of the analysis of the govern- 
ing equations. The values are also important because they form the basis 
for predicting the changes in the properties in those cases where the 
plasma was seeded with tungsten. 

Electrical Conductivity 


The mobility, p , of a gas is defined in terms of the drift velocity 
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of the electrons, u and the applied electric field, E 

“ - § ( 42 ) 

and the electrical conductivity, a, in terms of the current density, J, 
is given by 

J = aE = — (43) 

M- 

If it is assumed that an electron loses all of its momentum in each col- 
lision, then 

qE - m e uv e (44) 

where q is the elementary electric charge, is the mass of an electron, 
and is the electron collision frequency. From equations (43) and (44) 
the average drift velocity is 

u = — (45) 

m v ' 

e e 


and the current density is 

J = qN u (46) 

where N is the number density of electrons. Substituting equation (45) 
into equation (46) 
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q*N E 


m V 

e e 


(47) 


or from equation (43) 




m v 
e e 


(48) 


The electron collision frequency can be expressed in terms of (the 

effective momentum transfer collision cross section between electrons and 
species i of gas j), the velocity of the electrons v g , and the number 
density of species i of gas j N^. 


v . - f Pi v e 


(49) 


i — i i 

where the mean quantity v Q . may be approximated by v Q . when Q , is 

e ei e ex. ex 


taken at v . Then 
e 


v 

e 


v Q J - 

e ] I x ex 


(50) 


If the electron velocity distribution is Maxwellian 


v 

e 



(51) 


where k is the Boltzmann constant and is the electron temperature. 


Then 



The Saha equation predicts the number density of electrons from the number 


density of neutrals of gas j , N. 


0 


3, 

>4 


N J = 


i e * (. P . (si^))’ 


“o< T > 


(53) 


where h is Planck's constant, cp^ is the ionization potential of gas j, 

AE is the lowering of the ionization potential due to increased charged 
particle interactions found in the high temperature high density environ- 
ment of a plasma, and u^ and are the partition functions for singly 

ionized and neutral atoms defined by 

* 





-E /kT 
s' 

g s e 


(54) 


s=l 


u; 



S s 


-E /kT 
s' 
e 


where g and g + are the statistical weights of the neutral and singly 
s s 

ionized atoms at excitation level s (s = 1 ground state) and E is the 

s 

bound energy of level s . Substituting equation (53) into equation (52) 


and simplifying 
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2uJ(T) 




a = 


h 2 ZShJ qL 
j i 1 


For the case of argon only, equation (55) becomes 


2 / ra e kTlT ey rI Ar b \ I 2u + (T) ( ~ 


°Ar = 


* 5 Nf £ 


(55) 


(56) 


Figure 4 is a compilation of results from the theoretical works of De- 

97 90 90 30 31 

voto, Yos, Cann, and Cambel, and the experimental work of Emmons, 

30 33 

Lin et al., and Lau. The empirical values used in this study are 
determined by an analytical expression that best represents the available 
experimental and theoretical data. 

For the case of argon seeded with submicron tungsten particles, 
the summation over j in equation (55) must now include both argon and 
tungsten. Performing this summation and factoring out the expression for 
pure argon conductivity one obtains 


a Ar+w . °Ar 


1 + 


,,w w Ar 

K o u + u 0 

\| Ar w Ar 

N N 0 u o u + 


exp 




1 + 


S N w Q w 
i i ei 

S N Ar Q Ar 

i 1 


( 57 ) 


Air w 34- 

Substituting in values for <p and cp from Drawin and Felenkok and 



Electrical Conductivity 



Figure 4. Theoretical and Experimental Values for Electrical 
Conductivity of Argon at High Temperatures 
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w 

Ar Ar 35 U + 

average values for u + and from Olsen and taking — 

36 u 0 

Laun equation (57) becomes 


= 1 from 


Ar+w 


1 + 


N X r N 

-~-J exp [7.78/kT] 




i + 


E N? Q w . 
I i ^ei 


(58) 


The electron collision cross sections for tungsten have not been 

calculated or experimentally measured, but based on the cross sections of 

37 38 

other metals and on the predicted values of uranium and similar metals, 

it was assumed that (£ Q^)/(£ N^ r Q^5) « 1. This assumption is sus- 

pect for tungsten number densities greater than one percent of the argon 
number density. 

The effect of the seed can be seen in Figure 5 which shows the 
electrical conductivity of seeded argon from equation (58) compared with 
the conductivity of pure argon. 

Thermal Conductivity 

A compilation of argon thermal conductivity data is plotted in Fig- 
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ure 6 showing the theoretical work by Devoto, Yos, Admur and Mason, 
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Cann, and the experimental work of Emmons, and Knopp and Cambel. 

An analytical expression that yields the empirical values used in 

the numerical analysis was fit to the experimental data. It is felt that 

the experimental values are more reliable since they are easy to measure 

and are in good agreement, while the theoretical basis for thermal con- 
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ductivity is quite complex. 

27 

Devoto conveniently treated the contribution of the electrons to 
the thermal conductivity of the plasma separately, so the effect of the 
addition of tungsten seed was assumed to be an increase of the electron 
contribution to the thermal conductivity, due to an increase in the elec- 
tron number density. The effect of the tungsten collisions on the heavy 
atom thermal conductivity was neglected because of the low tungsten par- 
tial pressures. For evaluations of the case of tungsten seeded argon, 
the thermal conductivity due to the heavy atoms was assumed the same as 
that of pure argon. The contribution of the electrons to the thermal con- 
ductivity was calculated by taking the value of the electron component 
of the thermal conductivity of pure argon at the temperature required to 
produce the same number density of electrons as was present in the seeded 
case. Figure 7 is a plot of the thermal conductivity assumed in the analy 
sis and a plot of the thermal conductivity of the tungsten seeded argon 
for tungsten number densities used in the experiment where the electron 
contribution is added at a teirq>erature that yields an equivalent electron 
number density for the case of pure argon. This approach to accounting 
for the tungsten is considered to be as good an approach as possible at 
the present time due to lack of data on heavy metal collision cross sec- 
tions . 

Radiation Source Strength 

The values used in the numerical analysis for the radiation source 
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strength term were taken from the d.c. arc data of Emmons modified below 
9000°K to bring them more in line with the theoretical estimates of Hom^ 
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and the more recent induction discharge data obtained by Stokes. It 

was felt that the induction arc data were more reliable at the lower tem- 
42 

peratures. Stokes points out that "An important difference between the 

induction discharge experiments reported here and the arcs used for most 

of Emmons work is that stable, much larger plasmas can be obtained at 

low peak temperatures with induction heating." Pridmore-Brown,^ in his 

numerical analysis, did much the same things at the suggestion of Eckert 

who felt that Emmons’ values were too high at low temperatures. Miller 
21 

and Ayen arrived at a radiation source strength term for their two- 
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dimensional analysis by extrapolating the data of Evans and Tankin to 
lower temperatures which yielded values significantly below those of 
Emmons at temperatures below 10,000°K. 

31 42 

Figure 8 is a plot of the experimental data of Emmons, Stokes, 

43 44 

Evans and Tankin, and Schreiber and Hunter along with the theoretical 

values derived by Horn,^ and the values predicted by the Kramer-Unsold^ 

theory. The empirical values chosen for this work are also plotted. The 

radiation source strength term was not modified for the addition of tung- 
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sten seed. Work by United Aircraft on tungsten seeded plasmas with 
similar operating characteristics shows less than a 10% increase in power 
radiated for tungsten partial pressures on the order of 0.005 atm. This 
assumption may not be valid for higher tungsten partial pressures. 

Effect of Magnetic Field 

The magnetic field strength in the coil for the case of no plasma 
represents a maximum, since the field strength within the plasma is re- 
duced. Assuming an rms current of 100 amps in the work coil and an effec- 
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tive turn per cm of 0.52 (verified by experiment — see Appendix A), the 

magnetic field intensity would be 52 amp- turns /cm, which corresponds to a 

-3 2 

magnetic induction of about 5 x 10 webers/m , or 50 gauss. The work of 
27 

Devoto has shown that it would take a magnetic induction of one to two 
orders of magnitude greater to cause any change in the transport proper- 
ties. Further, the results of the numerical analysis and the experimental 
47 48 

work of Eckert and Trekhov et al. show the magnetic induction inside 
the plasma to be an order of magnitude lower than the field strength at 
the plasma radius. 
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CHAPTER III 

INSTRUMENTATION AND EQUIPMENT 
RF Plasma Generator 

The radio-frequency plasma generator and associated gas flow system 

was designed and built by Dr. Robert A. Benns as part of his dissertation 

research. A detailed description of the equipment and theory of operation 

49 

are found in his dissertation. 

The generator operates at 4 MHz and utilizes an air cooled Eimac 
4CX-35000C tetrode with an output power of up to 85 kW. The torch and 
operating controls are located in a copper lined room with a viewing port 
in the wall that allows spectroscopic observation of the plasma from the 
adjacent diagnostics room. The copper shielding is necessary to insure 
that stray rf fields are not present in the diagnostics equipment. The 
plasma generator operating parameters are listed in Table 1. 

The plasma is contained in a 40.64 cm long quartz tube with an 
inside diameter of 5.08 cm and a wall thickness of 1.75 mm. A l/4 inch 
water cooled copper work coil surrounds the torch. The work coil is ap- 
proximately 9 cm long with a coil pitch of approximately 2/3 turn per cm. 
The pitch of the work coil is not constant but has been varied to allow 
a large enough space between the second and third turns from the bottom 
to permit scanning of the entire diameter of the plasma without coil 
interference. The exact geometry of the coil and torch is drawn to scale 
in Figure 33 of Appendix A. 
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Table 1. Plasma Generator Operating Parameters 



Operating Value 

Maximum Rating 

dc plate voltage 

16,500 volts 

20,000 volts 

dc grid voltage 

- 400 volts 


dc screen voltage 

750 volts 

2,500 volts 

dc plate current 

6.2 amperes 

15 amperes 

dc grid current 

.14 amperes 


dc screen current 

.50 amperes 


peak rf plate voltage 

15,900 volts 


peak rf plate current 

9.3 amperes 


peak rf grid voltage 

500 volts 


peak rf grid current 

.28 amperes 


input impedance 

1,785 ohms 


output impedance 

1,892 ohms 


total plate input power 

102,000 watts 


plate dissipation 

28,000 watts 

35,000 watts 

grid dissipation 

14 watts 

500 watts 

screen dissipation 

200 watts 

1,750 watts 

grid driving power 

70 watts 


plate output power 

74,000 watts 


plate efficiency 

72.5 percent 




35 


A water cooled copper electrode is inserted through a size 12 
rubber stopper in the bottom of the torch. The tip of the electrode is 
coated with tungsten. The electrode serves as ground for the arc welder 
during the starting process and also helps stabilize the vortex flow in 
the tube. Figure 9 shows the quartz torch with the stopper and electrode 
in place. The center of the electrode is open to allow a regulated amount 
of gas flow through the center of the tube. The primary vortex gas flow 
is from the two copper tubes inserted near the perimeter of the stopper. 

The ends of the copper tubes are closed off and a small hole drilled in 
each provides flow tangential to the tube wall inclined 45 degrees in the 
upward direction. 

The gas flow system is designed to allow any desired combination 
of pure argon and aerosol to enter either the tangential flow copper tubes 
or the center flow through the water cooled electrode. In practice, very 
little flow could be sustained through the center electrode without plasma 
extinction, thus only a small pure argon flow was used which raised the 
visible plasma flame off the surface of the electrode. Samples of the 
gas were taken when tungsten aerosol was injected in the plasma to deter- 
mine the aerosol density. This was accomplished by operating the solenoids 
that allowed the vacuum tank to draw in a measured volume of aerosol 
through previously weighed filter papers. Two samples were taken for 
every aerosol run to insure a consistent aerosol and to give some basis 
for ascertaining the accuracy of this measurement. 

The submicron tungsten particles are injected into the plasma in an 
argon aerosol using the same gas flow system. This is accomplished by 
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allowing a set fraction of the pure argon to flow into a large test tube 
filled with dry tungsten powder. The injected argon agitates the reser- 
voir of tungsten powder on the bottom of the tube and creates an aerosol 
in the upper volume of the tube. The tube is mounted vertically on a vi- 
brator and the aerosol is drawn off the top of the tube and mixed with 
the pure argon that bypassed the tube. The major drawback is that it 
takes approximately two minutes for the aerosol density to stabilize. 
Because the plasma coupling has such a large dependence on the aerosol 
density, stabilization of the screen and grid currents of the power ampli- 
fier tube is a good indication of aerosol density stabilization. Figure 
10 is a drawing of the test tube configuration. The tube is mounted ver- 
tically on a vibrator. 

Diagnostics Instrumentation 

The plasma diagnostics instrumentation is contained in a room 
adjacent to the plasma generator with a viewing port in the wall between 
the two rooms. A schematic diagram of the instrumentation and optics 
used for spectroscopic observation of the plasma is given in Figure 11. 

The optics are mounted on a four by six foot Ealing optical table. 
Magnetic optical rails and carriers are used so that the optical compo- 
nents can be easily adjusted. A lens with a focal length of 743 mm 
focuses an enlarged image of the plasma on the entrance slit of the spec- 
trometer. The scanning mirror, mounted on a custom made precision scanner, 
rotates without vibration. This causes the plasma image to move across 
the entrance slit of the spectrometer. A potentiometer is mounted with 
amplification gears on the scanning mechanism and the variable resistance 
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is calibrated to provide a spatial position trace on the visicorder out- 
put. Spatial filters are used to insure that no stray reflected light 
from the plasma is reflected through the viewing port. A laser is used 
to locate the optical path and allow the spatial calibration of the scan- 
ning mechanism. This is accomplished by noting the visicorder indication 
from the scanning potentiometer when the laser is at the torch radius. 

The spectrometer is a McPherson Model 2051 scanning monochromator 
with a one meter focal length and an f 8.7 plane diffraction grating. 

The grating is 102 mm x 102 mm with 1200 lines/mm and is blazed for 5000 A. 
An EMI 9558QA photomultiplier tube is mounted on the exit slit. The en- 
trance slit width and the exit slit to the photomultiplier are set at 
50 p for all continuum observations. The light entering the spectrometer 
falls on only a small part of the grating and remains entirely on the 
grating during the spatial scan. The photomultiplier has 1000 V applied 
which is in the linear range of the tube. 

The output of the photomultiplier tube and the scanning potentiom- 
eter is monitored by one of three picoammeters and recorded by a visi- 
corder which uses light sensitive paper and pinpoints of light to leave 
a permanent trace of the picoammeter output. The response time of the 
system is limited by the visicorder which has a full scale response time 
of approximately .02 second. 

The absolute calibration of the system is accomplished with a 
tungsten strip lamp and regulated power supply traceable to the National 
Bureau of Standards. The lamp is mounted as shown in Figure 11 and a re- 
motely movable mirror, identical to the mirror in the plasma room, is used 
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to provide a substitute optical path. Half of a quartz tube is placed in 
the path of the lamp to simulate the torch. Figure 12 is a plot of the 
NBS tungsten strip lamp at 35 amps vs. wavelength and Figure 13 is a plot 
of the absolute sensitivity of the optical system vs. wavelength. 

The coil current is measured with a Tektronix P6021 high frequency 
current probe and Type 134 current probe amplifier fed into a Tektronix 
Model 597 cathode ray oscilloscope. The current probe saturates at 1 amp 
rms so a Tektronix Model CT-5 high frequency step-down transformer was 
mounted around the plasma work coil and was set to a 1000:1 step-down 
ratio. Tektronix calibrated the transformer and current probe to an accu- 
racy of three percent. A Polaroid camera attachment was used on the 
cathode ray oscilloscope to record the current trace. 



Lamp Intensity (watts/cm SR) 


2 






Sensitivity (p,amp cm SR/watt) X 10 
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Figure 13. Sensitivity of Plasma Diagnostics System 
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CHAPTER IV 

EXPERIMENTAL PROCEDURE 

A typical experimental run consists of from six to ten spatial 
scans of the plasma at different power settings. The first five scans 
are generally pure argon. The tungsten generator is turned on and the 
flow system is adjusted to inject tungsten aerosol into the argon for the 
remaining scans. The flow is allowed to stabilize and the tungsten flow 
is not adjusted for the scans of plasma containing tungsten. 

A new quartz torch is used for each run since the high power levels 
required for sustaining a plasma with tungsten aerosol usually destroys 
or permanently discolors the tube used in the previous run. 

The torch assembly is mounted inside the work coil and the gas flow 
lines are attached. Hie entire assembly is then checked for gas flow 
leaks since any air entering the torch will make plasma ignition virtually 
impossible. New filter papers are dried out, weighed, and inserted in 
the sample holders. 

The new torch is spatially calibrated by setting the scanning 
monochromator to the helium-neon laser wavelength (6328 A) and opening the 
exit slit. The optical path inside the diagnostics room is checked to 
insure the laser beam is centered on the mirrors and lens and the path 
it traces during the scan is centered in the rectangular slits of the 
spatial filters. The position of the laser beam inside the quartz torch 
is calibrated to the visicorder output of the rotating potentiometer that 
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is gear mounted to the scanning mechanism and has a small d.c. voltage 
applied to it. 

The diagnostics room equipment is warmed up and the plasma genera- 
tor start-up procedure is initiated. The plasma start-up check list is 
given in Appendix B and takes approximately forty minutes to complete. 

The scanning monochromator is set to the desired continuum wavelength, 

4315 A for most runs, and the slit widths checked at 50 p, . The high 
voltage to the photomultiplier is increased to 1000 volts. 

The plasma is ignited with a d.c. power supply operating at 40 
amperes. The plasma generator can only be raised to approximately 4 kW 
plate discharge power before ignition because the tank circuit is tuned 
to the impedance of the argon plasma. The rf driver voltage is raised 
until the tube screen current is near the value for automatic over-current 
trip-off. An arc is drawn from the welding rod to the grounded water 
cooled center electrode. The generator power is raised after ignition, 
and the center electrode is then ungrounded. 

The plasma is run for approximately five minutes at a plate dis- 
charge power level of approximately 15 kW to allow the amplifier tube to 
warm up and the quartz torch wall temperature to stabilize. 

The plasma is then adjusted to the desired power level and the 
d.c. plate voltage and current are recorded. The desired picoammeter 
scale is determined and recorded while observing the center of the plasma. 
The visicorder is then turned on and the plasma scanned. The current in 
the work coil is observed on the cathode ray oscilloscope during the scan 
and the current amplifier is set to scale and recorded. A Polaroid photo- 
graph is taken of the sine wave with the scale factor and the run and scan 
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numbers recorded on the back of the photograph. 

The power levels required for the seeded plasma are higher than 
for the pure argon plasma so the power must be raised during the switch- 
over to a seeded plasma. The vacuum pump is turned on and a vacuum drawn 
on the sample tank. The data scans of the seeded plasma are accomplished 
in the same manner as those of the pure argon. The low power scans are 
done first and the power raised for each run. At high power levels the 
tube will discolor. The visicorder plot of each scan is checked and 
asymmetry of the scan is used as an indication of tube discoloration and 
thus ends the run. 

Two aerosol samples are taken with the plasma running. The vacuum 
of the sample tank is recorded and the solenoid valve that opens the tank 
to the flow system is opened for approximately 10 seconds. The valve is 
then closed and the tank vacuum recorded. A sample is then taken through 
the remaining filter paper with the solenoid valve open for approximately 
twenty seconds. 

The filter papers are removed and weighed and the aerosol density 
is calculated from the weight change of the filter papers and the volume 
of argon drawn into the sample tank. 

The visicorder output is cut up into the individual scans and each 
scan is checked for symmetry. The run and scan number are recorded on 
the visicorder output. The photographs of the current of the work coil 
are measured and, based on the amplifier scale, the peak to peak and root 
mean square currents are calculated and recorded. 

The visicorder output with the recorded picoammeter scales serve 
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as the raw data for the spatially observed intensity profile used in the 
data analysis section that follows. The calibration of the observed in- 
tensity is accomplished with a NBS standard tungsten strip lamp as de- 
scribed in Chapter III. A typical visicorder trace is shown in Chapter 


V. 
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CHAPTER V 
DATA ANALYSIS 
Introduction 

There is a wide range of diagnostic methods available to the 

experimenter for measuring the electron temperature of a plasma. These 

methods fall chiefly into three categories; probes, and the active and 

passive uses of electromagnetic radiation. Probes, such as those de- 

50 

veloped by Tonks and Langmuir, are historically significant but not 
practical for high density plasmas due to deterioration of the probes 
and perturbations the probe may introduce. Thompson scattering of pho- 
tons by electrons yields information on electron temperature and density 
and ion temperature. The cross section for Thompson scattering is so 
small that lasers with powers of hundreds of megawatts are required to 
overcome the background due to the bremsstrahlung radiation. Passive 
methods that involve observation of the emitted radiation include tech- 
niques utilizing x-rays, microwaves, and radiation in the visible, infra- 
red, and ultraviolet range. Relative intensities of spectral lines, line 
shapes, or the intensity of continuum radiation may be employed. Meald 
and Wharton' 5 * give a good general tabulation of useful diagnostic tech- 
niques and the temperature and density ranges over which the methods may 
be utilized. 

The electron temperatures (on the order of 6000-9000°K) and number 
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densities (10 21 -10 22 m~ 3 ) of the argon plasma under consideration favor 
spectroscopic techniques in the visible range. Equipment requirements 
are minimized in the visible range and the emission intensity is high. 

The temperature dependence of the absolute continuum emission was chosen 
as the basis for electron temperature measurement in favor of several 
other available techniques. The chief advantage of this method is the 
high degree of dependence the continuum emission coefficient has on tern - 
perature in this range, thus minimizing the resulting experimental error 
in temperature determination. This technique does require an absolute 
calibration of the system. 

Methods involving the relative intensities of spectral lines require 
only relative calibrations of the optical system but errors are introduced 
in the form of experimentally determined atomic transition probabilities 
and possible self absorption of lines. These transition probabilities are 
readily available, 52,53,5 ^ but the estimated error in some argon lines is 
as high as 30%. This technique also requires that two or more spectral 
lines at different wavelengths be observed simultaneously while spatially 
scanning the plasma. A Boltzmann atomic plot of the natural log of the 
relative intensities of many lines versus the upper level electron energy 
associated with the transition that leads to each line reduces some of the 
errors caused by the uncertainties in transition probabilities. The slope 
of a straight line through the data points is a function of temperature. 

The accuracy of this method increases with additional lines and with lines 
spread over a wide range of electron upper level energies. Unfortunately, 
the slope of the line is not highly temperature dependent and, for argon 
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lines in the visible range a wide range in upper level energies requires 
a wide range in the observed spectrum. A Boltzmann plot was made from 
observations at the center of the plasma using lines between 4000 A and 
8000 A with poor results. 

The measurement of absolute line intensities also introduces experi- 
mental problems. Either the line shape must be determined at many spatial 
points with a high resolution spectrometer, or a slit width that introduces 
the entire line radiation must be determined and the continuum background 
taken into account. Self absorption at the line center may also be a 
problem. 

The continuum measurements were taken at the 4315 A wavelength and 
checked by comparison with three runs at 5000 A. A photoelectric trace 
of the argon spectrum at an average temperature of 7500°K from 4200 A to 
4400 A with a tungsten seed density of 18 x 10 ^ gm/cc is given in Figure 
14 and shows the continuum at 4315 A to be 15 A from any argon or tungsten 
spectral lines. 

Analysis Equations 

There are two independent but experimentally inseparable contribu- 
tions to the continuum radiation emitted from a plasma in local thermo- 
dynamic equilibrium (LTE). Bremsstrahlung (radiation from free-free 
transitions) is the result of the Coulomb interaction of unbound electrons 
with positive ions. The kinetic energy of the electrons is transformed 
into radiation as the electrons are accelerated by the electric and mag- 
netic fields present in the plasma. Free-bound (recombination) radiation 
is due to the transition of free electrons into bound states; the reverse 
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of photoionization . 

Kramers 55 first calculated the emission coefficient due to Brems- 
strahlung for a Maxwellian electron velocity distribution and a plasma 
containing only ions of equal charge. The emission coefficient expressed 
as energy radiated per unit frequency is: 


ff 16 tt e 6 z 2 N e N. G ff 
Gv 3c 2 <6TTmk)fc 


{-M 


(59) 


with the quantities defined as follows: 
e = electronic charge 
z = level of ionization 

= number density of electrons 
= number density of singly ionized atoms 
G^ = free-free Gaunt factor 
c = speed of light 
m = rest mass of electron 
k = Boltzmann constant 
T = electron temperature 
h = Planck’s constant 
v = frequency 

The Gaunt factor corrects the semi -classical photoionization cross 
section for hydrogen like systems to the exact values calculated through 
wave mechanics. The Gaunt factors, averaged over the Maxwellian electron 
velocity distribution, have been calculated by Karzas and Latter. 5 ^ The 
values of the free-free Gaunt factors are of order 1 and are usually taken 
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to be unity at the lower end of the visible spectrum, since the free-free 
contribution to the continuum is small compared to the free-bound continuum* 
The free-bound emission coefficient for a Maxwellian electron ve- 
locity distribution involves a summation of emission coefficients over all 
possible energy states to which a free electron may undergo transition. 
Unsold^ demonstrated that the summation could be replaced by an integral 
for frequencies that involved transitions to sufficiently dense levels. 

Unsold defines this cut-off frequency v in terms of the energy difference 

& 

between the ionization potential and the bound state for which higher 
levels are spaced closely enough to allow integration. The emission co- 
efficient for a one-electron system at frequencies below v is: 

s 
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Biberraan et al."* ’ corrected the emission coefficient for other hydrogen- 
like gases by introducing a correction factor C(T,v) or 
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i {l - exp (- £*)} §(T,v) (61) 


where g z ^ is the statistical weight of the ground state of the singly 

ionized atom and u is the partition function of the singly ionized atom 
z 

defined by: 


i K 

= 1 **,, 
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Bauder 60 used a cut-off frequency of v = 9.5 x 10 14 sec' 1 , corres- 

s 

ponding to a wavelength of \ = 3160 A in his work on high density d.c, 

61 62 

plasma arcs. Schluter * has calculated § factors for argon which are 
nearly temperature independent below 14,000°K. 

The total continuum emission coefficient is then the sum of e fb and 

v 

ff 

% ♦ Taking the free-free Gaunt factor to be unity and expressing the 
emission coefficient as energy radiated per unit wavelength by multiplying 
equations (59) and (61) by c/A^ we have: 
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(63) 


where 

a _ hv _ hc_ 
kT AkT 

Evaluating the constants in equation (63) and using the data from 
34 

Drawin and Felenbok for the statistical weight (g = 6) and partition 

Z )1 

function (u^(T) = 5,519) for argon over the temperature range of interest , 

equation (63) becomes (for c, in — - atts , ) ; 

A _ 4 

ster. cm 


A > 3160 A 


— 1.65 x 10 


-19 


e i 6g (A ,T) - a , -<y~| 

^2 ^ L 5.519 (1 * e > + e J (64) 


where 

N g = electron number density, cm" J 

= singly ionized number density, cm' 3 
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T — electron temperature, °K 
X — continuum wavelength, k 

The number density of the electrons is equal to that of the singly 
ionized atoms below approximately 20,000°K and can be calculated from the 
Saha relation, equation (53) in Chapter II. 

The neutral particle number density can be approximated from the 
ideal gas law for slightly ionized argon by 

P - N R T (65) 

o 

where R is the universal gas constant and P is the pressure which is very 
nearly 1 atmosphere pressure in the plasma. 

The number densities for argon neutrals, singly ionized atoms, and 
electrons at one atmosphere are given in Figure 15 as a function of tem- 
perature as calculated from equations (53) and (65) • The continuum emis- 
sion coefficient as a function of temperature given by equation (64) for 
argon at 4300 k and one atmosphere is plotted in Figure 16* 

The assumption of local thermodynamic equilibrium (LTE) requires 
that the population densities of the electrons be determined exclusively 
by particle collisions and that all species have a Maxwellian velocity dis- 
tribution at the same temperature* For bound levels, the distributions are 
given by the Boltzmann and Saha equations. Griem presents a criterion 
for LTE based on a collisional rate which is ten times larger than the 
radiative transition rate. If we apply this validity criterion estimate to 

the experimental conditions present in this work, we find that electron 

22-3 

densities of about 5 x 10 m are required which makes the LTE assumption 
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6 ^ 

suspect below about 11,000°K* Shumaker and Popenoe used a wall stabil- 
ized arc for studies of equilibrium in argon and concluded that departures 

22 -3 

from LTE occur at electron densities below 5 x 10 m Their data in- 
cluded no correction for self-absorption of the lines studied* Olsen^ 
used an argon arc burning freely between cooled metal electrodes to com- 
pare temperatures determined from measured emission coefficients of three 
species of plasma radiation and concluded that his atmospheric pressure 
plasma was in LTE at a temperature of 10,000°K. He states that "Discrep- 
ancies between temperatures determined from atom and ion line intensities 
have been shown to have been caused by self -absorption of the atom line 
in the outer layers of the plasma*" 

Work on rf argon plasmas by Scholz and Anderson^ indicates that an 
rf plasma may be in LTE at much lower number densities than a d*c* arc 

because of the larger plasma radii and lower axial and radial temperature 

64 

gradients present* Shumaker and Popenoe discussed the relevance of 
plasma diameter and temperature gradients in their report on d.c. arc 
studies* 

By comparing an arc with a black body in which complete LTE 
exists one can recognize three essential differences which may 
cause departures from equilibrium* First, in the arc much of 
the radiation escapes so that LTE will be achieved only if radi- 
ative depopulation rates are negligible compared to rates of 
collisional population and depopulation. Second, in the arc 
steep temperature and number density gradients exist which may 
generate diffusional fluxes of such magnitude that collisional 
(and radiative) processes are unable to maintain equilibrium 
distributions* And third, since electrons absorb energy from an 
electric field much more rapidly than ions do the kinetic tem- 
perature of the electron gas will be higher than that of the 
heavy particle gas which obtains its energy then mainly by col- 
lisions with the electrons* Since collisional processes will 
be dominated by electron collisions it turns out that all of 
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these effects should decrease rapidly in importance with 
increasing electron density . . «♦ The importance of the 
second effect should also depend upon arc channel diameter 
and radial position in the arc.^ 

Scholz and Anderson"* investigated LTE in an rf argon plasma operat- 
ing at 4.2 MHZ at pressures of 0.01, 0.1, 0.5, and 1*0 atmosphere. They 
used four emission lines and the continuum centered at three different 
wavelengths to show that the central core of the rf argon plasma was in 
LTE at 1.0 and 0.5 atm and very nearly so at 0.1 atm. They state that 
this 

. . . equilibrium is due in part to the existence of low (i.e. 

| VT | ^ 530°K/cm) axial and radial temperature gradients. At 
these pressures the plasma is in an apparent state of LTE and 
the thermal limit extends down to the ground state. The LTE 
electron number density criterion given by equation (1) [Griem] 
appears to be conservative for this plasma. “ 

The problem of relating the radial emission coefficient to the ob- 
served intensity profile is simply a geometric one for optically thin 
plasmas and the standard Abel inversion technique may be used. The prob- 
lem is more difficult when reabsorption inside the plasma cannot be ne- 
glected. Griem^ suggested a self-consistent iterative procedure for ob- 
taining emission coefficients in the case of moderate optical depths. 

The absorption coefficient is defined by: 

< 66 > 

where B(\ ,T) is the blackbody radiation. 

The plasma is assumed axially symmetric and divided into N zones 
of radii R/N , 2 R/n, . . ., and R as indicated in Figure 17. The tempera- 
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ture, emission coefficient, and absorption coefficient are considered 
constant in a given zone. The reduction in intensity of radiation I 
after passing through a zone of constant K(\,T)and distance x is given 
by: 

1 = 1 e- Kx (67) 

o 

If that zone is also emitting radiation with the same emission coefficient 
throughout that zone, the observed intensity will be: 

_ -Kx P* -K(x-£) i C / zq \ 

I = I e +J e e Vb/ d5 (68) 

o 

or 

I = I e" Kx + £ [l - e~ Kx ] (69) 

O K 

The intensity of the radiation observed at a given spatial position, 
y^, as shown in Figure 17, is the sum of the radiation emitted from each 
zone minus the radiation attenuated by the zones between each emitting 
zone and the surface* 

A computer program was written to relate the emission coefficient, 
absorption coefficient, and temperature of each zone to the observed in- 
tensity of the radiation emitted from the surface of the plasma. 

Computer Program — MAD 

The computer program, MAD (Modified Abel Diagnostic code), utilizes 
the spatially observed surface radiation intensity as raw data fed in at N 
discrete points according to Figure 17. The program calculates the radi- 
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ally dependent temperature, emission coefficient, and absorption coeffi- 
cient for a plasma that is not optically thin. A typical visicorder 
trace, as shown in Figure 18, is reduced to give values of the observed 
radiation in watts/sec m based on the output current of the photomulti- 
plier tube and the calibrated sensitivity of the system. The radius of 
the plasma is found from the radiation intensity trace and the spatial 
position trace. The radius is assumed to be that position where the in- 
tensity drops to zero. 

The data fed in are first curve fitted with a third order spline 
fit subroutine and interpolated to give observed intensity values for N 
zones. The program then calculates values of x(n,i), which is the thick- 
ness of zone i when observed at chordal position y , shown in Figure 17. 
This distance is 

x(n,i) = ^[(l - " l 1 " n) ] ' Z < 7 °) 

j=i+l 

During the first iteration, the emission coefficients are calcu- 
lated for each zone from the observed intensity I(y^) a nd the values of 
x(n,i), assuming the plasma to be optically thin, i.e*, the absorption 
coefficient K is set equal to zero for all zones. Starting from zone one 
and working toward the center, the emission coefficient of each zone is 
calculated using the emission coefficients of previously calculated zones 
from: 

n 

I Cy n ) =2 ^ e(i) x(n,i) (71) 

i=l 
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or 


:(n) 


n-1 

= [ I( V ' 2 1 e(1) x(n > 1) ] 2^7 

i=l 


(72) 


where e(i) for i < n have been calculated previously. 

The temperature and number density for each zone can be calculated 
from the emission coefficient and equations (53), (64), and (65). Values 

for the number density are interpolated from data given by Olsen's GASCOMP 

5 /4 b /T o 

program using an N = aT ' exp ' dependence between values 1000 K apart. 

The a and b coefficients are calculated for each 1000°K interval and the 
error is less than one percent over the temperature range of 4000-14, 000°K. 
The correct temperature and number density for a given emission coeffici- 
ent are then arrived at by binary search iterative procedure. The first 
approximation of the absorption coefficient is calculated by dividing the 
emission coefficient by the blackbody radiation at the temperature of the 
given zone as indicated by equation (66) . 

After the first iteration, a new set of emission coefficients is 
calculated based on the observed intensity data and the first approxima- 
tion to the absorption coefficients. Equation (71) is modified to include 
the absorption coefficients: 


i-1 


f( y n ) = Y {i<(i) ^ 1 ■ exp (~ x ( n > i ^ { exp ( Y _k ^ x ( n >j)) ( y3 > 

j=l 


1=1 


n 


+ exp ( Y x ( n >j) + Y -k ^ x ( n >j))}} 

j=l j=i+l 
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This calculation is best accomplished on the computer by defining 
a set of optical arrays based on equations (67) and (69). Let e(i)A.(n,i) 
represent the emission from a zone of length x(n,i) and constant emission 
coefficient e(i) and absorption coefficient K(i). From equation (69) 

A(n ,i) = [l - exp (-K (i) x(n,i))] (74) 

Defining B(n,i) as: 

B(n,i) -K (i) x(n,i) (75) 

then exp B(n,i) represents the attenuation of radiant energy passing 
through a zone of length x(n,i) according to equation (67). 

Arrays C(n,i) and D(n,i) are defined in terms of sums of B(n,i) 
such that exp C(n,i) represents the total attenuating power of all of the 
zones between the zone emitting radiation e(i)A(n,i) and the surface for 
zones on the near side of the plasma, and exp D(n,i) represents the equiv- 
alent for zones on the far side of the plasma: 

i-1 

C(n,i) = £ B(n , j) (76) 

j=l 

n n 

D(n,i) = Y B (n,j) + £ B(n,j) 

j=l j=i+l 


( 77 ) 
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Now the radiation emitted from zone i and along chordal position 

n is 

I = e (i)A(n ,i) ( 78 > 

The radiation from this zone reaching the surface is then 

I = e (i)A(n ,i) exp (C(n,i)) (79) 

for zones on the near side of the plasma and 

I = e (i)A(n ,i) exp (D(n,i)) (80) 

for zones on the far side of the plasma. Equation (73) then reduces to: 
n 

I(y n ) = Y, e(i)A(n,i) [exp (C(n,i) + exp (D(n,i))] (81) 

i=l 

The emission coefficients can then be calculated starting from zone 1 and 

proceeding inward by: 

n-1 

I(y ) - ; e(i)A(n,i) [exp (C(n,i)) + exp (D(n,i))] 

n 

e (n) = (82) 

A(n,n) [exp (C(n,n)) + exp (D(n,n))] 

where e(i) for i < n have been calculated previously. 

These new emission coefficients are used to calculate new tempera- 
tures and absorption coefficients. The iteration terminates when the 
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changes in emission coefficients, absorption coefficients, and temperatures 
for all zones are less than one percent of the values in the previous 
iteration* 

This analysis provides a self-consistent radial distribution of 
temperature, emission coefficient, and absorption coefficient for the dis- 
tribution of the observed intensity for moderately self-absorbing plasmas. 

A listing of the program is given in Appendix C. 

Error Analysis 

Chordal Average Temperature 

The experimental errors associated with the temperature measurements 
are small compared with the experimental magnetic field intensity measure- 
ments because of the high temperature dependence of the emission coeffi- 
cient. If we define a chordal average temperature as the equally weighted 
sum of the temperature of each zone divided by the number of zones, and the 
chordal average emission coefficient as the value of the emission coeffi- 
cient determined by that temperature, the chordal average emission coef- 
ficient multiplied by the plasma diameter will yield the value of the 
emitted intensity, I(y c ) , along a chord through the center of the plasma. 

It is estimated that the photoelectric trace of the observed inten- 
sity can be read within 10% and the sensitivity of the system as shown in 
Figure 13 is known within 15% giving a total uncertainty in the observed 
intensity of 25%. The absorption coefficient is neglected in estimating 
the error in measured chordal average temperatures since the theory and 
experiment show the plasma is optically thin. The uncertainty in observed 
intensity results in an uncertainty in the average emission coefficient 
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of the same value since 


so 


and 


e = Ky c )/R 


(83) 


« = lr “c = “ c /R 

C 


(84) 



? Ke 



(85) 


The uncertainty in the chordal average temperature can be esti- 
mated from the observed intensity error by 


ST 

AT 5z 
T T 


Ae 


g 

T I 


( 86 ) 


The values of e and — as functions of temperature can be obtained 
from equation (64) or estimated from Figure 16. Over the temperature 
range of 7000°K £ T £ 9500 d K, a 25% uncertainty in observed intensity re- 
sults in a chordal average temperature uncertainty of less than 2%. The 
uncertainty in experimentally measured chordal average temperatures is 
much less than that of the experimentally measured magnetic field inten- 
sity. Thus, only the error bars associated with the magnetic field inten- 
sity measurements are shown in the results section. 

Magnetic Field Intensity 

The Tektronix high frequency current probe, amplifier, and step 
down transformer are calibrated to within 3%. It is estimated that the 
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magnetic field intensity can be calculated from the coil current to within 
an accuracy of 2% based on the experimental measurements reported in Ap- 
pendix A. The current recorded from the photographs of the trace on the 
oscilloscope is estimated to be correct within 10%, giving a total uncer- 
tainty in the reported magnetic field intensity of 154. 

Tungsten Aerosol Density 

The tungsten aerosol density was measured twice during each experi- 
mental run and the uncertainty in the reported density is estimated to be 
20%, based on the average discrepancies in measured values. 
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CHAPTER VI 
RESULTS 
Introduction 

A total of 59 experimental runs resulting in temperature profiles 
were made. Approximately 50 more runs were made that did not yield re- 
sults due to equipment failure. The failure was most often in the plasma 
generating equipment. 

The first 27 runs were made without measurement of the coil current 
in order to develop a degree of confidence in the plasma generator and the 
diagnostics instrumentation. Three runs were made using the continuum at 
5000 1 and at 4315 JL to ensure that the resulting temperature profiles 
were not dependent on the continuum wavelength chosen. The plasma spec- 
trum was recorded from 3500-8000 A and analyzed to ensure that no unknown 
contaminants were in the plasma. Twelve of these first 27 runs included 

various concentrations of tungsten aerosols ranging in number density 

3 -2 

ratios N /n, from 2.5 x 10" to 3.7 x 10 to ensure that the effect of 
w' Ar 

tungsten could be measured and to gain experience in operating the plasma 
generator with a seeded plasma. High speed, 16 mm motion pictures of the 
plasma were made at 1000 frames/sec to check for instabilities in the 
plasma that might not be observable to the naked eye. A small pulsation 
of the plasma at a frequency of 60 Hz could be detected which caused a 
variation in the plasma radius of about 107«. This resulted from the 60 Hz 
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ripple in the coil current from the rectifier circuit. 

There were 32 runs made with temperature profiles recorded at 
various power levels. The coil current was measured along with the plate 

dissipation power. Eleven of these runs were made with tungsten to argon 

-3 -3 -2 

number density ratios of 4*3 x 10 , 6.1 x 10 , and 3.6 x 10 . 

Exper imen tal Results 

Argon 

The observed intensity for four representative runs at magnetic 
field intensities of 33.1, 46.0, 55.2, and 77.2 amp-turns/cm are illus- 
trated in Figure 19. Figure 20 is a plot of the resulting radial distri- 
butions of the calculated emission coefficients. The absorption coeffi- 
cients were also calculated in the analysis and ranged from about 10 6 

1 _5 _1 

cm for the low power case to 2 x 10 cm for the highest power case, 

indicating that the plasma was optically thin. 

The temperature profiles that result from the radial emission co- 
efficients are presented in Figure 21. While the observed intensity 
ranged over an order of magnitude for the various power levels, the re- 
sulting temperature profiles varied only about 10% from the lowest power 
case to the highest. This is due to the very high temperature dependence 
of the emission coefficient. 

The four temperature profiles presented in Figure 21 are compared 
with the theoretical analyses in Figure 22. The solid line indicates 
the experimental results; the dotted and dashed lines represent the nu- 
merical results described in Chapter II. The plasma radius was reduced 
in one case to the experimentally observed radius and the thermal con- 
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Figure 19. Observed Intensity Profiles for Four Representative 
Magnetic Field Intensities 
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Figure 20. Calculated Emission Coefficient Profiles for Four 
Representative Magnetic Field Intensities 
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Figure 21, Experimental Temperature Profiles for Four Representative 
Magnetic Field Intensities 
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ductivity increased as explained later. The off-center peak temperature 
of the numerical analysis is much more pronounced than the results of the 
experiment. The off-center maximum of the temperature profile is a result 
of the skin heating by the electromagnetic field and radiation cooling of 
the plasma. When radiation is neglected, the numerical analysis yields a 
temperature profile with the maximum temperature occurring at the center 
of the plasma. The experimental results also have an off-center tempera- 
ture peak but the difference between the center temperature and the peak 
temperature is much less than predicted by theory. This is in keeping 
with the results of other experimenters as shown in Chapter I. The most 
probable reason for this is the complex turbulent flow pattern present 
in the experimental plasma that is not accounted for in the theoretical 
analysis. This turbulent mixing flow will tend to flatten out the tem- 
perature profile. Neglecting axial convection in the numerical analyses 
may also account for the difference since the axial temperature gradient 
will be highest at the radial position where the peak temperature occurs, 
tending to flatten the temperature profile. The magnetic field intensity 
is increasing at the axial point, where the plasma was observed, as shown 
in Figure 33 of Appendix A, and the peak temperature increases with in- 
creasing magnetic field intensity. An approximate analysis was made with 
an assumed axial temperature gradient calculated from the derivative of the 
applied magnetic field intensity and the change in the radially dependent 
temperature with applied magnetic field intensity as predicted by the nu- 
merical analysis. This analysis shows that the correction term is an order 
of magnitude lower than the radiation term. This is not in agreement with 
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21 

the theoretical predictions of Miller and Ayen. Their two-dimensional 
analysis of temperature profiles in a torch with an assumed flat axial 
flow shows a definite flattening of the temperature profiles at the begin- 
ning of the energy addition region and a gradual increasing of the off- 
center maximum temperature at increasing axial positions. 

The effect of the turbulence is to increase the effective thermal 
conductivity. The numerical analysis shows a marked decrease in magnetic 
field intensity near the tube wall while it appears that the vortex flow 
decreases the radius of the plasma, thus delaying the drop in magnetic 
field values. The numerical analysis was run with the thermal conductivity 
values doubled and the radius reduced to match the experimentally observed 
radius. The results are shown in Figure 22. The profile shapes are in 
better agreement with the experimentally observed profiles, indicating 
that the turbulent flow may indeed cause a significant increase in the 
effective thermal conductivity due to mass convection. 

Additional insights can be gained if the chordal average temperatures 
are plotted as functions of applied magnetic field strengths as shown in 
Figure 23. The chordal average of the numerical temperature profiles re- 
mains in good agreement with the experimental chordal average throughout 
most of the range of magnetic field intensities. The uncertainty in the 
numerical solution due to the uncertainty in the electrical conductivity 
is reflected by the cross hatching. 

The analytical approximation of Roninskii et al. predicts tempera- 
tures several thousand degrees higher than does the numerical analysis. 

This is due to the radiation term which is neglected. At lower tempera- 
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tures, radiation plays a less important role, and the two techniques con- 
verge. If the radiation term is set equal to zero in the numerical tech- 
nique, the two solutions are nearly identical. In this case the maximum 
temperature occurs in the center for both techniques. 

Figure 24 is a plot of the total linear power input and the power 
radiated as a function of magnetic field intensity predicted by the nu- 
merical analysis. The radiated power represents 50% of the total power 
input for a magnetic field strength of 75 amp-turns/cm and drops to 35% 
for a field strength of 30 am-turns/cm. 

The numerical analysis predicts that a plasma can be sustained at 
magnetic field strengths as low as 15 amp-turns/cm, but experimentally the 
limit was determined to be about 30 amp-turns/cm. This can be seen in Fig- 
ure 23, where the experimental points tend to diverge from theory at about 
30 amp-turns/cm. This can also be attributed to flow. The plasma is 
literally "blown out" at low power levels by the flow. 

Argon Seeded with Tungsten 

The argon plasma was seeded with tungsten at number density ratios 
-3 -3 -2 

of 4.3 x 10 , 6.1 x 10 , and 3.6 x 10 . The greatest effect on the 
temperature profiles was with an N^/N^ of 6.1 x 10 . This is most likely 
because the collision cross section of tungsten can no longer be neglected 
at the higher number densities of tungsten. Temperature profiles for 

t 

three different magnetic field strengths and a tungsten number density 
ratio of 6.1 x 10 are shown in Figure 25. The temperatures are reduced 
approximately 1000°K below the pure argon case. Figure 26 is a plot of 
these same three profiles compared with the results of the numerical 
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Figure 24. Numerical Analysis Predictions of Total Linear Power 
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analysis. The transport properties were altered in the numerical analysis 
as discussed in Chapter II. Figure 27 is a plot of the chordal average 
temperatures of the plasma as a function of magnetic field strength for 
the various tungsten partial pressures. The temperatures decrease with 
the addition of seed as predicted by the theoretical analysis with gener- 
ally good agreement until the seed reaches an N /N. of 3.6 x 10“ 2 . The 

W A.I" 

experimental temperatures begin to rise for higher seed concentrations 
while the temperature calculated from the theory continues to decrease. 
This is because the theory neglects electron collisions with tungsten in 
the calculation of electrical conductivity. 

Parametric Study of Numerical Analysis 
The uncertainties in the transport properties of both seeded and 
unseeded argon and the experimental parameters that can be varied make it 
worthwhile to examine the effect the data input to the numerical analysis 
has on the predicted tempera ture profiles. Figures 28, 29, 30, 31, and 
32 each shows the temperature profiles at four different magnetic field 
intensities with the electrical conductivity, thermal conductivity, fre— ■ 
quency, radius, and radial mass flow rate varied over the indicated ranges. 
The effect of the electrical conductivity can be seen in Figure 28. 
The values were varied by translating the empirical curve with respect to 
temperature by 1000°K to either side of the values normally used in the 
analysis. This effectively bounds the available data. The electrical 
conductivity terms enter the governing equations through both the energy 
balance equation and the equation for the induced electric field. The 
induced electric field is lowered for higher values of electrical conduc- 
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tivtty. This reduces the power input term for the energy balance equa- 
tion, oE , since the increase in o is more than offset by the decrease 

w 

in E * The reduction in the power input term causes a uniform reduction 
0 

in temperature, though the temperature values do not seem to be very 
sensitive to variations in the electrical conductivity values. 

The temperature profiles are not very sensitive to the values 
chosen for the thermal conductivity either. Figure 29 shows the results 
of varying the values of K by 25%. An increase in thermal conductivity 
causes the temperature profile to flatten as might be expected, but does 
not lower the average temperature in spite of the fact that more energy 
is transferred through the torch wall. This is most likely because the 
reduction in the maximum hot spot temperature lowers the total radiated 
energy more than the temperature increase in the central region increases 
it. 

The frequency of the applied electromagnetic field effects its pene- 
trating ability. At higher frequencies, the power input is over a nar- 
rower band near the tube wall, resulting in a lower center region temper- 
ature but a higher temperature near the wall. This effect can be seen 
in Figure 30. 

The results of variations in the plasma radius are not so easy to 
predict. An increase in radius requires a larger net power input to keep 
the same values of the applied magnetic field intensity at the plasma 
radius. The surface area for heat transfer to the wall also increases. 

The hot spot temperature remains approximately the same as seen in Figure 
31, though the central temperatures drop with increasing radius since the 
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induced electric field decreeses as the penetration distance increases* 

A study of the effect of a radial mass flow was made since there 

22 

is some experimental evidence that indicates such a flow may exist* 
The results are shown in Figure 32* The hot spot temperature is not 
affected in magnitude or position. The center region temperatures are 
increased with an inward flow toward a sink in the center. 

The most surprising result of these studies is the very low sen- 
sitivity of the temperature profiles to any of the inputs > including 
the magnetic field intensity. 
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CHAPTER VII 

CONCLUSIONS AND RECOMMENDATIONS 
Conclusions 

The one -dimensional system of equations made up of the energy 
balance equation. Maxwell's equations, and Ohm's law are adequate to 
quantitatively describe the chordal average temperature inside an rf 
plasma over a wide range of magnetic flux intensities but do not describe 
the profile shape as well as might be desired. This is a possible result 
of the complex turbulent flow patterns inside the torch which are not ac- 
counted for in the analysis. There does not appear to be any reasonable 
solution at present to this problem, though experimental and theoretical 
temperature profiles would most certainly be in better agreement for ex- 
perimental plasmas designed to minimize this effect. An effort was made 
to simulate the effects of turbulence by increasing the thermal conduc- 
tivity values and using the experimentally observed radius in place of 
the torch radius in the numerical analysis. Agreement with experimentally 
observed temperature profiles was improved. 

The neglect of the radiation term in the energy balance equation 
is not justified. The results of the numerical treatment and experiment 
show that the temperatures encountered at the magnetic flux intensities 
necessary to sustain a plasma are too high to neglect the effect of radia- 
tion. The closed- form analytical approximations that neglect radiation 
do yield interesting insights into the problem, most notably with regard 
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to the degree of coupling between the electrical and thermal conductivi” 
ties and the resulting temperature predictions, but do not accurately 
predict the max imum temperatures encountered or the temperature profiles. 

The numerical solutions are not highly sensitive to any of the 
experimentally controllable inputs or to the electrical or thermal con- 
ductivity values chosen. This lack of sensitivity was experimentally ob- 
served with respect to the magnetic flux intensity and conductivity values 
but equipment limitations prevented variations of the torch radius or rf 
driving frequency. The addition of tungsten seed significantly increased 
the electrical and thermal conductivities of the argon but lowered the 
temperatures of the plasma, as predicted by theory, by only about 1000°K 
for corresponding magnetic field intensities. 

Recommendations 

An experimental arrangement that minimizes the effect of turbulent 
flow would form a better basis for comparison with theory. Extending the 
experimental measurements to two dimensions would be a valuable addition 
to those numerical techniques which attempt to take into account the 
axial variation of the energy addition region. Utilizing gases with 
marked differences in transport properties would help establish the 
ability of the theory to handle the problem. A more productive effort 
might be made in the direction of simplifying the calculations necessary 
to determine the chordal average temperatures resulting from the many 
variables open to the experimenter. A modification of the existing closed 
form solutions that takes into account radiation would be a valuable ad- 
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dition. A detailed comparison with existing numerical solutions over 
a wide range of cases would lend credence to the results of this modified 
closed form solution. 
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APPENDIX A 

COIL AND MAGNETIC FLUX GEOMETRY 

The magnetic field intensity inside the work coil is a function of 
the work coil geometry and is proportional to the current in the work 
coil. Figure 33 is a full scale drawing of a cross section of the work 
coil and torch with the plane of observation noted. 

The magnetic field intensity at any point z along the axis can be 
calculated from the work coil current, I, by 


H 


2 

i 


a 2 I 


2 (a 2 + x i )^ 


(87) 


where a is the work coil radius and x i is the axial distance from point z 
to each coil position, as noted in Figure 33. 

Defining 


2 

Y(z) = 2 y 

i 2 (a + x £ ) 


( 88 ) 


Then 


H(z) = I Y(z) 


(89) 


The calculated values of y( z ) are shown in Figure 33. These values 
were experimentally verified by passing a known d.c. current through the 
coil and measuring the magnetic field intensity inside the coil with a 
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Bell Model 640 Gaussmeter. Figure 34 is a plot of experimentally meas- 
ured points taken at the observation plane used for the plasma tempera- 
ture profiles and a line corresponding to the theoretically calculated 
value of Y = .52 at that point. The agreement appears excellent. 






APPENDIX B 


PLASMA START-UP CHECK LIST 


Pre-Start 

Pad Input Switch 
Arc Welder 
Main Water 

Capacitor Drain Switch 

New Filter Papers 

Spectrometer 

Torch Alignment 

Argon Supply 

Bypass Flow 

Torch Leaks 

Visicorder Spatial Alignment 
Spectrometer 
Spatial Scan 
Ground Switch 

Start 

Capacitor Drain Switch 
Pad Input Switch 
Arc Welding Rod 
Coil Cooling Water 


off 

off 

off 

down 

check 

6328 k 

check 

check 

clean 

check 

record 

4315 1 

centered 

on 

up 

on 

clear 

on 


Interlock System 


all off (left) 



100 


Pad Input Starter 

on 

Rectifier Switches 

on (right to left) 

Plate Voltage Delay 

on 

Screen Voltage Control 

full c*c.w. 

Grid Voltage Control 

full C.C,W. 

Grid and Screen Power Supply Circuit 
Breakers 

on 

RF Driver 

full c.c.w. 

Filament Voltage 

full c.c.w. 

lube Blower Circuit Breaker 

on 

Aux. Blowers and d.c. Arc Circuit Breakers 

on 

Filament Circuit Breaker 

on 


1 kw Amp. Circuit Breaker 
RF Driver — four switches 


Grid Voltage 
Screen Voltage 



Filament Voltage 
Plate Voltage 
Plate High Voltage 
Plate Voltage 
Grid Voltage 



Screen Voltage 
Filament Voltage 
Grid Voltage 
Screen Voltage 
Filament Voltage 


on 

1 4 

on 2 3 

200 V 

200 V 

2 V 

lower until zero light on 
on 

2500 V 
300 V 
300 V 
3.5 V 
400 V 
400 V 
5 V 
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("Grid Voltage 

5 min 

Screen Voltage 
J?i lament Voltage 
Cascade Interlock Lights 
Ready Light 
^Filament Voltage 
Switching Network 
Turnoff 
Grid 
Plate 
Screen 
RF Driver 


Plate Voltage 
RF Driver 


500 V 
700 V 
7 V 

check on 
check on 
9.7 V 


on (right) 
on (right) 
on (right) 
on (right) 

c.w. until screen current noted, 
then fyll c.c.w. 

5000 V 

c.w. until screen current noted, 
then full c.c.w. 


Plate Voltage 
RF Driver 

Argon Pressure 
Argon Bypass Flow Rate 
Torch Blowers 
Arc Welder 
RF Driver 
Ignition 


7500 V 

c.w. until screen current noted, 
then full c.c.w. 

40 psi 

3 " 


on 


on 


c.w. until screen current 150 ma 


Draw Arc with Welding Rod to 
Water Cooled Center Electrode 
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Immediately after Ignition 
Arc Welder 

RF Driver c.w. until 

Ground Switch 
Duct in Place 
Warm up Torch Slowly 
Normal Shut Down 

Simultaneously 

Switching Network Turnoff 
RF Driver 
Immediately 

Switching Network Screen 
Plate 
Grid 

Filament 

Plate Voltage 
Filament 
Screen Voltage 
Grid Voltage 
Plate High Voltage 
Delay 

Rectifier Switches 


off 

screen current 100 ma 
off 


off (left) 
full c.c.w. 

off (left) 
off (left) 
off (left) 
off (left) 

2000 V 
0 V 
0 V 
0 V 

off 

off 

off (left to right) 

stop 

off 


Pad Input Starter 
Argon Flow 

Grid and Screen Supply Circuit Breakers 


off 



Filament Circuit Breaker 
RF Driver--four switches 
1 kW Amplifier Circuit Breaker 
Torch Blowers 
Water 
^ 15 Min 

Aux. Blowers and d.c. Arc Circuit Blowers 

Tube Blower Circuit Breaker 

Water 

Pad Input Switch 
Main Water 
Capacitor Drain 
Unintentional Plasma Shut Down 
S imu 1 taneous ly 

Switching Network Turnoff 
RF Driver 
Inmediately 

Switching Network Screen 
Plate 
Grid 

Filament 

Re-Start 

Plate Voltage 
Overcurrent Trips 


off 

-- 4 1 

off 3 2 

off 

off 

1/2 flow 


off 

off 

off 

off 

off 

off 

off (left) 
full c.c.w. 


off (left) 
off (left) 
off (left) 
off (left) 

3000 V 
reset 


Filament 


9.7 V 
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Ground Switch 
Switching Network: 

Turnoff 

Grid 

Plate 

Screen 

RF Driver 

Plate Voltage 
RF Driver 

Plate Voltage 
RF Driver 

Arc Welder 
RF Driver 


on 


on (right) 
on (right) 
on (right) 
on (right) 

c.w, until screen current noted, 
then full c.c.w. 

5000 V 

c.w. until screen current noted, 
then full c.c.w. 

7500 V 

c.w. until screen current noted, 
then full c.c.w. 

oh 

c.w, until screen current 150 ma 


Follow Normal Ignition Procedure 








APPENDIX C 


MODIFIED ABEL INVERSION COMPUTER PROGRAM 



c 

c 

c 

c 

#-> 

c 


INTEGER CoDEtY,Yl 
INTEGER COUNT 

real ic, kappa 
real KOLD 


real NE 

dimension 

dimension 

dimension 

dimension 

dimension 

DIMENSION 


IC(5o) »X(50tSl),EC(sO) , KAPPA (50) »3(50) 

0lC(26),RE(l)*AC(26)rBC(26),CC( 2 6),DC(26)»CCC(2S) 
A(50 r 51) »RX(S0»5l) »C(50r51) |OC50»51> 

TOLD (50) #KOLO(50 > * ECO (50) 

XUU) »YUCi) »YNR'(26> rOlCN(26) 


FIRST CARD CONTAINS RUN CODE » PLASMA RADUjS(CM) »NUMBER OF ZONES 
USED IN ANALYSIS (BETWEEN 4 AND 50 ) f WAVELENGTH (ANGSTROMS ) t NUMBER 
DATA POINTS FED IN (BETWEEN 4 *ND 25) 
formated (i3,fio.2»ho»fi o.iriio) 

IF THIS CARD IS BLANK PROGRAM ASSUMES END OF DATA SETS 
99 READ (5*000) CODEf R»N» WL»ND p 
N =40 


i 


1 


1 


1 

i 


\ 


OF 

1 

i 


IF(CODE.Eq.O) go to 920 
WRITE (S*3o 1)CODE'R> Nr WL»NDP 
WRITE ( i f Boo ) CODE» R»N» WL»ND P 

c ' NEXT READ IN OBSERVED INTENSITY (WATTS/CM**3-S p > FOR EACH ZONE»THE 
c OUTSIDE ZONE IS Y=1 . ONE DATA p OlNT OM EACH CARD WITH THE ZONE NUMBER 

c and the: intensity ASd3f ireip. 3) 
c the program assumes the intensity is 
C ZERO AT THE RADIUS of THE PLASMA 


DO 50 Y= 1 ,N 
TOLD(Y)=0,0 
KOLD(Y)=0.0 
50 CONTINUE 
KNOT=C 


COUNTSO 


WRlIc(6f B02) 

DO 100 K=1»NDP 
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READ ( 5 f 8C 3 ) Y»OlC(Y) 

WRlTE(6»a04)Y»0lC(Y) ! 

WRITE 1 1 » 816) Y# OIC ( Y ) , 

100 continue 

DO 101 I=i, MDP 
Cl = I-.l 

C2=N0P S 

YNR(I)=R*(C1/C2) ■ 

OICNd)=OIC(NDP-I + l) 

101 continue ! 

YNr ( NDP+1 ) =R 

" YNr(NDP+2)=R*(J+(1/C2) ) 

01CN(UDP+1)=0.0 ! 

0ICN(NDP+2) =0.0 

C^uL COEFFf YNR,OICN,NDP+2tl*0»2f 0 ♦ RE » P Q » AC » 3C » CC , DO CCC » Xu» YU) , 
WRI7E16, 8i3) 

DO 102 Y = 1 » N 

C1=N . - . 

C2 = Y 

C3=R*( (C1-C2)/(CJJ ) "... . 

--- -CALL ?NTRP4(C3'»C4»1) 

IC(Y)=C4 

WRlTE(6f6l4) Y , C3» IC ( Y ) 

WRITE ( l r 614 ) Y»C3»IC(Y) 

102 CONTINUE 

c CALCULATE GEOMETRIC, matrix ARRAY _ _ . ... .1 

- X t N, N+l) =0 » 0 

DO 105 K=Ni»l»-l 
DO 105 I=K..l»-l 
X(KrK+l)=o.O 
Z2=0.0 

IF(I.E0,K) GO TO 104 
11 = 1+1 

DO 105 J=I1,K 
Z2=Z2 + XC< f J) 
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1C3 CONTINUE 

104 YY=I-1 
1Z-K 

ZYZ=N „ 

X(«.r 1 ) = (R)*( (1-YY/ZYZ) **2-(l-ZZ/ZYZ)**2>**0.5-Z2 

105 CONTINUE 

c calculate Emission coefficients -assume plasma optically thin 
EC(i)=(lC(l))/(2*X(i»l ) ) 

IF(EC(1).LT.0.0> EC<1>=0,1 

•• uo 250 K=2.N 

SUM=0.0 

Ji=K-i 

DO 210 J=lrJl 
SUM=S’JM+X ( K > J ) *EC ( J) *2 , 0 

210 continue 

EC(K>=(IC(K)-SUrt>/<2*X<K,K> > 

IF{EC(K) .LT.0.0) EC(K)=0.0 
250 CONTINUE 
c CALCULATE TEMPERATURE 
300 CONTINUE 

DO 350 K-Nt 1»-1 

• STH=10999,0 " " 

GTl=4001,o 

310 5T=(GTH+3TU/(2*0) 

NE=ARND(St) 

GEc=EC£q(GT»NE» WL) 

IF(GEC.ST.EC(K) ) GTH=3T 
IF{GEC,LE,EC(K) ) GTL=GT 
IF((GTH-STL).GT. 05,0)50 TO 310 
T(K>=<GTH+$TL)/(2.0> 

350 CONTINUE 
C CALCULATE OPACITY 
400 continue 

DO 425 K=Nfl»“l 
3 (K) =BLK3DY (VfLf T(K ) ) 
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t^wppfUK) =EC IK) /a(K> 

<+25 continue 
500 CONTINUE 
JJ =0 

J SO ■ ■ . 

COUNT-COUnT +1 
00 525 Y=1,N 

* PC£nT=(T( Y)-TOLO(Y) )/(T( Y) ) 

PCENO= ( kappa ( Y ) -KOLO < Y ) ) / ( KAPPA ( Y) ) 

pcene=<ec(y)-e:cO(Y) >/(ec(y) ) 

IF(PC£NT. SE. 0*005) J=1 
IF(PCEM0*SE. 0,005) JJ=1 
IF(PCENE,Gt. 0,005) J = 1 
TO L D(Y)=T(Y) 

KOLO( Y)=KAPPA(Y) 

ECo(Y)=EC(Y) 

525 CONTINUE 

IF ( J , EO , 0) GO TO 900 
IF(JJ.EQ.o) GO TO 900 

IF (COUNT, £0.20) 00" TO 900 ' " 

c calculate optical arrays 
600 continue 

K0ut=K0MT+1 
00.620 Y=1,N 
DO 620 I=\, Y 

A ( y, I i = ( 1-EXP ( -KAPPA (I)*X(Y»I) ) ) /( KAPpA ( I ) ) 
BX ( Y» I ) =-KAPPA ( I ) *X{ Y» I ) 

C(Y»I)=0.0 
OlY,U=0.0 
62n CONTINUE 

C ( 1 , 1 ) =0 • 0 

0(l,l)=BX(l,l) 

□0 650 Y=2»M 
DO 650 I = J,Y 
IF(I.EQ.I) go to 630 
SUm= 0.0 



11 = 1-1 

DO 623 L-l r X 1 
SUM=SU-UBx(Y»L) 

625 CONTINUE 
C ( Y» I) =SIN 
630 CONTINUE 
SU"=0.0 
DO 635 L=iVY 
SUM=SJH+Bx(YrL) 

635 CONTINUE 

IF { I *cQ, Y) GO TO 641 

11 = 1+1 ' 

DO 640 L.=n»Y 

~ SUM=SUM+3X(Y»L) ’ “ ~ ' - 

640 CONTINUE 

641 D<Y»I)=SUM 
650 CONTINUE 

C CALCULATE new EMMISSION COEFFICIENTS 
700 CONTINUE 

EC ( 1 ) = ( I c ( n j / ( A c rm +Ail, lt*EXP tux ( 1 V 1 > ) ) ~ “ 

IF(EC(1).LT.0.0> EC<1>=0.0 
DO 747 Y=2»N 

SUM=0 .0 

Y1=Y-1 

DO 737 L=1 t Y1 

SUM=SUM+(ECtL)*A(Y»L) )*(Ex p tO(Y»L) )+Ex p <C{ YrL) ) ) 

737 CONTINUE 

EC(Y>=(IC(Y)-SUM>/<A(Y»Y)*(EXP0(Y,Y) )+EXP<C(YrY>) ) ) 

IfTEC(Y) .LT.O'.O) EC(Y)=0,0 

747 CONTINUE 
GO TO 300 

800 FORMAT ( I3*F10*?»I t0»Fl0.1»H0) 

801 FORMAT ( lHl» 20X» f ANALYSIS PROGRAM FOR ARGOn PLASMA* »///» 

1 30 X, ’RUN CODE IS* ,!4»/r 

2 30X» ' PLASMA RADIUS IS* »F6V3* IX# *CM. * t /* * 


i 
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3 30X* ’THERE ARE»*I3*1X* »ZONES USED IN THE ANALYSIS**/* 

4 30X» ’ THE CONTINUUM IS OBSERVED AT ♦ »P7. 1 » IX* * ANGSTROMS’ * 

5 /, 30X, ’THERE WERE* * 14* IX* » DATA PolNTS’) 

802 FORMAT ( 1H1 » 20X» ’ OBSERVED jNTENSl TY« * /// » 1SX » ♦ ZONE’ » lOX* 

I ’INTENSITY (WATTS/SR-CM**4) ’*//) 

803 F0RMAr<l3,lPE19.3) 

804 rORMAT( 15X* 14* lOX* 1PE9.3) 

809 FORMAT<lHlf30X, ’CALCULATED' PLASMA PARAMETERS’*///* 

1 10X»’20NE’»5X* ’OBSERVED INTENSITY’ *5X* ’EMISSION COEFFICIENT’* 

2 5X* ’TEMPERATURE’ »5X* ’ABSORPTION COEFFICIENT* >5X* ’8LACK BODY RAD* 

3*/*l9X* » (WATTS/CM**3-SR> * *5X” <WATTS/CM**4-SR) ’* 

A" 5X* »OES. KELVIN » *5X»* ~ (CM-1) * ri4x* » ( wATTS/CM**3-sR) ♦ r//> 

310 FORMA TdOx* I3*5X>1PE14.3,7X*1PE17,3*9X’0PF9.1»10X* 

1 1PE12.3*10X,1PE12.3*/} 

813 FORMA f< 10X* * ZONE* ♦ lOX* ’ R* , 20X* ’ IC ♦ *///) 

814 FOrmAT(10x*I3*F14,3»F20*3) 

815 F0RMAT(I4,3F20,5) 

316'FORMAT{T3 f F10.1'> ~ 

817 FORMAT (120) 

900 continue 

c write out data 
WRITE(6,809) 

00 910 Y=N*i*-l 

WRITE (&» 810) Y , I c ( Y ) f EC ( Y ) * T C Y> » K APPA ( Y> * B < Y > 

WRITE ( 1 * 815 ) Y, EC(Y) »T(Y) *KAPPA(Y) 

910 CONTINUE 

WRXT£(6*617) KONT 
50 TO 99 
920 CONTINUE 

end * • ' - • 
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FUNCTION GAUNT (WL) 

C THIS FUNCTION CALCULATES the SAunT 
,c factor as a function of wavelength 
c (wl-amsstroms) emperically based on 

C- DATA BY DIETER SChLUTEr-Z . ASTROPHYSI K 
C VOL 6i PAGE 67-19&S-PANGE 3500-6500 A 

IF( (WL.LT.3500, > .OR. (WL.ST.&500. ) ) GO TO ? 

IF (WL.LT. 4800 , ) 00 TO 2 
IF( WL.LT.4900. ) GO TO 1 
GAunT=( 0. 5/600, >*( WL-4900. ) +1.5 

RETURN • ' — - 

1 GAijnT=(0. 4/100. )*(4900, -wL>+1.5 

return 

2 GAunT= ( 0. 75/1300. )*(WL-35oO. ) +1.15 
RETURN 

3 WRITE (6. 4) WL 

' "4 FORMAT ( lHlf IX r » WAVELENGTH OUT OF RANGE IN GAUNT»»F10, 2> 

GAyNT=2 . 0 

return - 

END \ 



FUNCTION BLKBOf(ArT) 

c this function calculates the blackbody 

C RADIATION IN(WATTS/SR-CM**3) for a 
C GIVEN WAVELENGTH (A-ANGSTROHS) AND A 
C GIVEN TEMPERATURE (T-DEG K). 

Cl=5.9Slti? 

C2si43,B6 
W-A/1000 . 0 
U=T/1000.0 

8LkbDY=( 1/(£XP{C2/U/W)-1) )*( (2*Cl)/(W**5) ) 
RETURN " ‘ 

end 
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FUNCTION ARND(T) 

C THIS FUNCTIONAL ROUTINE CALCULATES the 
_C ELECTRON and ION NUMPEr density IN CM - 3 

c (which are assumed equal) fop argon - ga'S” 

C AT ONE ATMOSPHERE AND WITHIN THE TEMPERATURE 
C RANGE OF 4000-14000 DEG K, THIS PROGRAM 

c interpolates between values that are iqoo deg k 

C APART ASSUMING AN N=A*T**(5/t|) *EXF(B/T) DEPENDENCE 

c the data points are taken from olsen# s gascomp 

C PROGRAM' REPORTED ON IN ARL 70-0048 MARCH' 1970 
DIMENSION DT(15> ,0N<15) 

IF( (T.LT. 4000. 0> .OR, (T.GE. 14000, 0) ) Go TO 201 

DO 100 1=1,15 
DT( I ) -I*10Q0, 0 

ON ( I ) = 0 , 0 

“100 CONTINUE • - 

ON ( 4 ) s 1.27E10 
0N{5) = 1.32E12 
ON ( 6 ) = 2.94E13 
0N(7) = 2.72E14 
DN( 8 > = 1,45E15 

- ON ( g ) r ' 5 . 34E 15"' * ' ' : 

ON(10)= 1.51E16 
ON ( 11 )= 3.46E16 
ON ( 12 ) = 6.62E16 
DN{13)= 1.08E17 
DN(14)= 1,50=17 

00 101 1=4,13 - • • — - 

IF{ (T.GE. DT(I) ) .AND. (T,LT,DT(I+ 1 ) ) ) GO TO 195 
' 101 CONTINUE 
GO TO 201 
195 A=DN(I) 

X=(ALOG(Dn(I+ 1 )/ON(I) ) >-(5*0/4.0>*{ALOG<DtU+1)/DT(I) ) ) 
B=(DT( I) *qT( I+l) *X) /(DT( I)“DT( 14-1 ) ) 

A=(DN(I) )/( (OT(I)**(5.0/4,0) )*(EXP(B/DT(I) ) ) ) 
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AR^D=A*T**(5.0/4.0)*EXO(8/T> 

return 

201 WRITE ( 6f 2o2) T 

202 FORMAT ( lHl , » TEMPERATURE 0(jT OF RANGE IN ArND T IS»»Flo.l> 
ARND=1.0El3 

return 

end 



FUNCTION £CEQ<T'N»WL) 

C THIS FUNCTION CALCULATES THE CONTINUUM 
C EMISSION IN WATTS/SR-Cm**4 OF ARSON 
C AS A FUNCTION OF tEMPERATURE(T-DEG.K) » 
C NUM3ER DENSITY OF IONS OR ELECTRONS 
C WHICH ARE ASSUMMED EQUAL (N-CM**-3> r 
C and WaVELEN3TH( WL-ANSSTROMS) , 

REAL KrNrNF 

H=if,14c-15 

ft—8*6l6E— r ) 

C-3« OtlQ 

V=C/WL 

U2=5.519 

X=EXP(-(H*V)/(K*T) ) 

Y=(6.0*3AUNT(WL) >/<U2) 

NF=N/X,0E5 

Z=(NF**2) /(SORT(T) *WL**2) 

ECEQ=U.65E-09*Z)*(Y*(1.0-X)+X> 

RETURN 

END 
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SU3R0UT I NE U5Q ( X t Y» NORSK , NY r 3 . CL , C » »« ) 
DIMENSION X(NY) »Y(NY) r3<NV> »CC(NY»MY) »C(NY) 
INTEGER P,Q,PQ 
PGrO 

m=norsr+i 

DO 10 0=1 , K 

3(®)=0. 

00 15 1=1, NY 

8tO)=B(0)+Y(I)*Xa)**(Q-l) 

15 CONTINUE 
DO 20 P=l,v 
CC(P*Q)=0. 

DO 25 1=1, NY 

CC tPrO)=CC (P,G)+X(I)** (P+0-2) 

25 CONTINUE 
20 CONTINUE 
10 CONTINUE 
DO 30 0=1, y 
DO 35 P=l,M 
L=a+(P-i)*v. 

C(U=CC{PrG) 

35 CONTINUE 
30 CONTINUE 

CALL SIMQ(C*B»M'KS> 

IF(KS.Eq.i) PQ=2 
RETURN 
—-END 
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SUBROUTINE: COEFFU* Y»NY»IFIT,NOROR,IPLOT,IR,rE,PQ» 
S A,b»C,D»CC»XU,YU) 

c 

DIMENSION X(NY) »Y(NY) 

DIMENSION A { NY ) » B ( NY ) » C l NY ) * D ( NY ) 

DIMENSION rC(NY'NY) 

DIMENSION XU (NY') »YU(NY) »RR(NY> 

INTEGER Pq 

data EPS/. 01/ 

GO TO 77 

ENTRY INTRP4(XI'YI»IFIT> | 

c : 

GO TO (31 r 32) » IFIT 

C ' i 

31 DO 20 1=1, NYP 
K=I 

IF(Xl.LT.X(I+l) ) GO TO 21 

20 CONTINUE 

21 Y I = A(K)*XI**3+B(K)*XI**2+C<K)*XI+D<K> 

“RETURN : • • • • — ' 

C 

32 Yl=0 • 

DO 34 J=l r MM j 

YI=YI+B( J)*XI**< j-l) 1 

34 CONTINUE ! 

return 

*** calculate fitting coefficients *** 

77 IF ( I PLOT . EQ • 5 .OR, IPLOT.eQ.8) GO TO 5^0 
IF(IPL0T.E0.7) GO TO 20C 

GO TO (4i,4R) t I FIT 

c v ! 

41 NYpzMY-l 


119 



DO 10 ,< = 1,NYP 

c compute constants 
X0=X(K-1> 

XK=x(K) 

X1=X<<+1) 

X2=X(<+2) 

YOzY ( K-l ) 

YK=y(X) 

Yl=YtX+l) 

Y2zYU+2) 

XC5=X0*X0 

XK$=XK*XK 

XlS=Xl*Xi 

X£s=X2*X2 

xoc=*os*xo 

XKc=XKS*XK 

X1c=X1S*X1 

X2c=X2S*X2 

XMq=XK-XO 

XMi=XK-Xl 

XM2=XK-X2' 

YMo=YK-YO 

YM i= YX-Yl 

YM2=YK-Y2 

xmos=x:<s-xcs 

XMlS=XXS-XlS 

XM2S=Xf<S-X2S 

XMoC=XKC-XOC 

XMlC=XKC-XlC 

XM2C=X’<C-X2C 
IF(K.E3.1) GO TO 11 
IF(K.EQ.NYP) go to 12 
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C FOR INTERNAL PANELS (K=2»NYP-1) 

A<«>= ( (YM0*XMi-YVIl*XMo)*(XMOS*XM2-XM2S*XM0) 

1 -< YM0*XM2-YM2*XM0>*(XMOS*XMl-XMls*XM0> > 

2 / r <XMOC*XMl-.XMlC*xMO)*(XMOS*XM2-XM2S*XMQ) 

3 -(XM0C*XV'2-XM2C*xM0)*(XM0S*XMl«XMlS*XM0) ) 

0<K)= ( ( YM0*XM2-YM2*XM0>-A(K)*(XM0C*XM2-XM2C*XM0) > 

1 / (XM05*X^2-XM2S*x SJ 10) 

CtK)= ( YMO-A^K)*XMOC-9(K)*XMOS > / IXMq) 

D(K)= YK - A(K)*XKC - B(K)»XKS - C(K>*XK 

- -- GO TO 10 " ' 

c for leftmost panel c k=i > 

11 A(K)=0, 

3<K>= ( YM1*XM2-YN2*XM1) / < XMlS*XM2-XM2S*XMI ) 

C<K)= ( YM1-3 ( K ) *XM1S) / ( XMl ) 

0(K)= YK-8(K)*XKS-C(K)*X< 

©0 TO 10 “ 

c for Rightmost panel (k=nyp> 

12 A ( K ) =0 . - 

B(K)= (YM0*XM1-YM1*XM0> / (XM0S*XM1-XM1S*XM0) 

C(K>= (YMO-B(K)*XMOS)/(“XmO) - — - - 

0(K)= YK-B (K ) *XKS-C ( K > *xK 

■ 10 CONTINUE • ~ - 

return 

42 MM=nOROR+i 

return 

c 

C LEftST SQUAREs-WITH 1/REL*ERR, WEIGHTING ' 

200 L=o 
MM = nOROr+i 

207 RMAX=0. 

DO 201 J=1,MY 

201 RMAX=AMAX1 (RMAx»RE( I) ) 

210 DO 202 1=1, NY 
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IF ( rEU) .LT.EPS) Rt(I)=EPS 
NUm=AMIN1 ( 10. »rMAX/RE ( I ) ♦EPS) 

IF(NUM.LE.O) So to 202 
00 203 

l=l+i 

--■XU(L)=X(I) 

YU(L)=Y('I) 

203 CONTINUE 
202 CONTINUE 

IF ( IPLOT . EO • 8 ) GO TO 520 
IF(L.LE.30) GO TO 205 

RMaxsRMAX/1,1* ~ 

WRITE16»299) L 

299 FORMAT ( 2lH WIDE WEIGHTING RANGE* I5» 14-H wTS REDUCED) 
L=0 

SO TO 210 

205 CA[.l LSQ(xU*YU,NORDR,UB,cC»CrPQ) 

— -RETURN • ■■■ — • 

c 

C LEAST SQUARES FIT TO HISTOGRAM DATA - 

C WITH OR WITHOUT l/REL.ERR. WEIGHTING 

500 MM = nORDr+i 
R=IR 

IF( IPL0T.EP.8) GO To 207 
520 LX=0 
NQr=0 

DO 1701 J=l ,MM ' 

Bt J)=0, 

DO 17U 1 = ) » NY 
NUM=I 

IF ( IPLOT, EO.B) NUM=AMIN1(iO, » RMAX/RE < I ) +EPS ) 
IF(NUM.LE.O) go to 1711 
DO 503 LMrl # NUM 

B(j)=3 (J)+( (r+i»>*(X(X+1)**(J+IR)-X(X)**(J+XR) )«Y<X) ) 

1 /(<FLOAT(U)+R)*lX(I+l)**(IR+l)-X(I)**(IR+n) ) 

* 
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503 CONTINUE 
17H CONTINUE 

DO 1702 Krl.MM 

c c ( 3 » fo = o . " 

DO 1712 1=1 » NY 
NUmsI 

IF(IPL0t.EQ. 6) NUM=AMINl(lO. »RMAX/RE<I>+EPS) 

IF(NU^,LE.0) GO to 1712 
DO 510 UMsl.NUM 

~ cc ( j> KTsCC ( J * KT* ~TT RTI ;y**21 (XTT*l'T**'t 0*1 RT-X CH **TJllR>1 

1 *(X(I+1>**U+IR)-X(I)**(K+IR)> ) 

z / t tploatiu)+R)* (Float (K' v+Rj^txt in) **(iRn> 

3 «xm**( IR+l> >**2 ) 

510 CONTINUE • 

1712 CONTINUE 

NBR=N3R1NUM — 

LX=J+ ( K-»l ) *MM 

: "Cl LX) =CC ( J, K ) •• ' ~ : T~~ 

1702 CONTINUE 

” 170r~CONTlNlJ£ ~ r •• - : 

IF (NBR.LE.225) GO To 11 li - 4^ - 

— rmax=r^a x/tt9 ~ t * i ^ A -“ 

WRITE (S, 299) NBR ■- - 

GO TO -520 — 

C 

'1111 CALL SIMQ(C»B»«Mp IT) " ' 

IF(IT.EQ.I) GO TO 999 _ . _ __ 

WRITE ( 6» 47) 

97 FORMAT (" 24K *!* SINGULAR MATRIX *** ) 

9999 RETURN 0 

999 RETURN - 

END 
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SUBROUTINE SIMQ(A» 9 #N»KS) 

- DIMENSION A ( 1 )> B(l) 

A - MATRIX OF COEFFICIENTS STORES T*E ROWS OF AN NXN MATRIX 
IN ONE COLUMN. THE N ELEMENTS 0 P THE FIRST ROW FORM THE 
FIRST N ELEMENTS OF THE A MATRIX8 jHE SECOND ROW FqRMS 
. - THe next n- ELEMENTS OF A ECT. ‘ ' ' ~~ 

A - MATRIX of coefficients stored .columnwise, these are 
destroyed in THE COMPUTATION, the SIZE of matrix a is 
N by N, 

' B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N). THESE AR£ 
REPLACED By FINAL SOLUTION VALUES, VECTOR X. 

H^ V'NUMBER 0'F~EQUATION^*AND VA-RTABLES.— - : -• 

KS - OUTPUT DIGIT, 

■ “ 0 FOR A NORMAL SOLUTION, * - - 

1 FOR A SINGULAR S^T OF EQUATIONS. 

TOljtO.O 

KS =0 

J -hi - - r -*• * ■ — - — 4 — • - - — 

DO 65 J=i,N ’ " ' ~ 

-UYSU+I" “ = 

JUsJJ+N+i 

BIG'A-O. ■ ~ • ' — - - 

IT=JJ-J 

T> 0 '* 30 ~r£ J , N" " — ~ ' ~ - - 

IU=IT+I 

IF ( A 3 SOIGA) - ABS(A(IJ) > > £ 0 , 30,20 ' 

20 BIGA=A(IJ) . .. 

IMax=I • • • 

30 CONTINUE .... 

IF ( 'A35 (STGA) - TOU 35, 35* 40 ' — - 

35 KS=l 

RETURN - - . ... - - - 

40 Ii=U+N*(J« 2 ) 

IT'slMAX-J- " ~ - - 



■ 00 50 K=J,N 

HrfU+M " 

ia=n+iT 

save=a(id -■ 

A<I1)=A{I2) 

A(I2>=SAVE 

50 A(I1)=A(U)/BIGA 

SAVE=8(IMAX5 

B(IMAX)=B( J) 

B(a)=SAVE/BISA 

IF(J-N)55,70»55 

- j-i) - — -r 

DO 65 I X = JY t N 

... ... tXdrWS+TjT"“ 

IT=J-IX 
80 60 UX=JY»N 
IXJX=N*( Jx-l)+lX 

• - J«JX=IX0X+IT 

60 A<IX<JX)=A<IXJX>-(A(IXJ)*A(JJX> ) 

— 6S~-fM I X') a 8i fxl -t © ( lOCd ) t~ ' 

70 NY=N-1 

fTsN*N- 

DO BO J=1»NY 

-IA = IT-J 

IBsN-J 

tesN . --- 

00 80 i<=l f J 

£H IBl ~8 ( Ig) -A 1 1 A ) *B ( IC ) 

J 80 lCsiC-1 * 

RETURN ■ 

END - 
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